Photinus pyralis (Linnaeus, 1767) is amongst the most widespread and abundant of all U.S. fireflies (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 non-feeding, 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).

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/elifesciences-publications/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

Download asset Open asset

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.

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 cytochrome-oxidase 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₂, 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 1—figure 2).

Appendix 1—figure 2

Download asset Open asset

1.3.2 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 dark-adapted 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 pre-cleaning 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 3

Download asset Open asset

Appendix 1—figure 4

Download asset Open asset

Appendix 1—figure 5

Download asset Open asset

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

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

1.5.1 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).

Appendix 1—figure 6

Download asset Open asset

1.5.2 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 µg 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 µL Qiagen QLE buffer (10 mM Tris-Cl, 0.1 mM EDTA, pH 8.5) yielding 17.1 µg of DNA (214 ng/µL) and glycogen (500 ng/µL). 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-gel-electrophoresis (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 µg 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 μg) 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

Download asset Open asset

Appendix 1—figure 8

Download asset Open asset

The P. pyralis genome assembly followed three stages: (1) a hybrid assembly using Illumina and PacBio reads, producing assembly Ppyr1.1 (Appendix 1.6.2), (2) Ppyr1.1 scaffolded using Hi-C data, producing assembly Ppyr1.2 (Appendix 1.6.3), and (3) Ppyr1.2 manually curation for proper X-chromosome assembly and removal of putative non-firefly sequences, producing Ppyr1.3 (1.6.4).

1.6.4 Ppyr1.3: Manual curation and taxonomic annotation filtering

1.6.4.1 Defining the X chromosome

Hi-C data was mapped and converted to the 'hic' file format with the juicer pipeline (v1.5.6) (Durand et al., 2016b), and then visualized using juicebox (v1.5.2) (Durand et al., 2016a). This visualization revealed a clear breakpoint in Hi-C linkage density on LG3 at ~22,220,000 bp. Mapping of Illumina short-insert and PacBio reads with Bowtie2 (v2.3.1) (Langmead and Salzberg, 2012) and SMRTPortal (v2.3.0.140893) with the ‘RS_Resequencing.1’ protocol, followed by visualization with Qualimap (v2.2.1) (Okonechnikov et al., 2016), revealed that the first section of LG3 (1–22,220,000 bp), here termed LG3a, was present at roughly half the coverage of LG3b (22,220,001–50,884,892 bp) in both the Illumina and PacBio libraries. 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 non-firefly 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.

Appendix 1—figure 9

Download asset Open asset

Appendix 1—figure 10

Download asset Open asset

Appendix 1—figure 11

Download asset Open asset

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) (Koren et al., 2017) 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.

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(Shen et al., 2016), 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).

Appendix 1—figure 12

Download asset Open asset

We annotated the coding gene structure of P. pyralis by integrating direct coding gene models produced from the de novo transcriptome (Appendix 1.9.2) and reference guided transcriptome (Appendix 1.9.3), with a lower weighted contribution of ab initio gene predictions, using the Evidence Modeler (EVM) algorithm (v1.1.1) (Haas et al., 2008). First, Augustus (v3.2.2) (Stanke et al., 2006) was trained against Ppyr1.2 with BUSCO (parameters: -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.0β) (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.

1.10.1 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 20-hydroxyecdysone 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).

Appendix 1—figure 13

Download asset Open asset

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

MethylC-seq libraries were prepared from HMW DNA prepared from four P. pyralis MMNJ males using a previously published protocol (Urich et al., 2015), and sequenced to ~36x expected depth on an Illumina NextSeq500. Methylation analysis was performed using methylpy (Schultz et al., 2015) Methylpy calls programs for read processing and aligning: (i) reads were trimmed of sequencing adapters using Cutadapt (Martin, 2011), (ii) processed reads were mapped to both a converted forward strand (cytosines to thymines) and converted reverse strand (guanines to adenines) using bowtie (flags: -S, -k 1, -m 1, --chunkmbs 3072, --best, --strata, -o 4, -e 80, -l 20, -n 0 [Langmead et al., 2009]), and (iii) PCR duplicates were removed using Picard (Picard Tools, 2017). In total, 49.4M reads were mapped corresponding to an actual sequencing depth of ~16x. A sodium bisulfite non-conversion rate of 0.17% was estimated from Lambda phage genomic DNA. Raw WGBS data can be found on the NCBI Gene Expression Omnibus (GSE107177). Previously published whole genome bisulfite sequencing (WGBS)/MethylC-seq libraries for Apis mellifera (Herb et al., 2012), Bombyx mori (Xiang et al., 2010), Nicrophorus vespilloides (Cunningham et al., 2015), and Zootermopsis nevadensis (Glastad et al., 2016) were downloaded from the Short Read Archive (SRA) using accessions SRR445803–4, SRR027157–9, SRR2017555, and SRR3139749, respectively. Libraries were subjected to identical methylation analysis as P. pyralis.

Weighted DNA methylation was calculated for CG sites by dividing the total number of aligned methylated reads by the total number of methylated plus un-methylated reads (Schultz et al., 2012). For genic metaplots, the gene body (start to stop codon), 1000 base pairs (bp) upstream, and 1000 bp downstream was divided into 20 windows proportional windows based on sequence length (bp). Weighted DNA methylation was calculated for each window and then plotted in R (v3.2.4) (R Development Core Team, 2013).

We synthesized a 5’ fluorescein-tagged (TTAGG)₅ 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).

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

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.

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 calibrated-fluorimetry 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 µL PBS. These tissues were chosen to avoid the ovary tissue. Once homogenized, 900 µL PBS, 1 µL Triton X-100 (Sigma-Aldrich), and 4 µL 100 mg/mL RNase A (QIAGEN) were added. The homogenate was incubated at 4˚C for 15 min, filtered with a 30 µm Cell Tries filter (Sysmex), and further diluted with 1 mL PBS. 20 µL 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 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), resulting in the final assembly (Alat1.3). The final assembly (Alat1.3) consists of 5388 scaffolds totaling 908.5 Gbp with an N50 length of 693.0 Kbp, corresponding to 96.6% of the predicted genome size of 940 Mbp based on flow cytometry (Appendix 2.4). Genome sequencing library statistics are available in Appendix 4—table 1.

Appendix 2—figure 1

Download asset Open asset

2.5.2 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 Mbp), representing 26.3% of the scaffold number and 1.3% of the nucleotides of Alat1.2, producing the final filtered assembly, dubbed Alat1.3.

Appendix 2—figure 2

Download asset Open asset

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.

cDNA libraries were generated from purified Total RNA (500 ng from each sample) using a TruSeq RNA Sample Preparation Kit v2 (Illumina) according to the manufacturer’s protocol (Low-Throughput Protocol), except that all reactions were carried out at half scale. The fragmentation of mRNA was performed for 4 min. The enrichment PCR was done using six cycles. A subset of nine libraries (BdM1, HeF1, HeM1, LtF1, LtM1, Egg1, Lrv1, PpEF, PpLM; Appendix 2—table 1) were multiplexed and sequenced in a single lane of Hiseq1500 101 × 101 bp paired-end reads. The remaining 23 libraries (BdM2, BdM3, HeF2, HeF3, HeM2, HeM3, LtF2, LtF3, LtM2, LtM3, WAF1, WAF2, WAF3, WAM1, WAM2, WAM3, Egg2, Lrv2, Lrv3, PpEM, PpLF, PpMF, PpMM) were multiplexed and sequenced in two lanes of Hiseq1500 66 bp single-end reads. Sequence quality was inspected with FastQC (Andrews, 2017).

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.

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

Ignelater luminosus is a member of the beetle family Elateridae (‘click beetles’), related to Lampyridae within the superfamily Elateroidea. The Elateridae includes about 10,000 species (Slipinski et al., 2011) (17 subfamilies) (Costa et al., 2010), 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; Costa et al., 2010). 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 (Costa et al., 2010), 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.

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.

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.

Appendix 3—figure 1

Download asset Open asset

The karyotype of male Puerto Rican I. luminosus (as Pyrophorus luminosus) was reported as 2n = 14A + X₁X₂Y (Virkki et al., 1984). The genome sizes of 5 male I. luminosus were determined by flow cytometry-mediated calibrated-fluorimetry of DNA content with propidium iodide stained nuclei by Dr. J. Spencer Johnston (Texas A&M University). The frozen head of each individual was placed into 1 mL of cold Galbraith buffer in a 1 mL Kontes Dounce Tissue Grinder along with the head of a female Drosophila virilis standard (1C = 328 Mbp). The nuclei from the sample and standard were released with 15 strokes of the ‘B’ (loose) pestle, filtered through 40 µm Nylon mesh, and stained with 25 mg/mL Propidium Iodide (PI). After a minimum of 30 min staining in the dark and cold, the average fluorescence channel number for the PI (red) fluorescence of the 2C (diploid) nuclei of the sample and standard were determined using a CytoFlex Flow Cytometer (Beckman-Coulter). The 1C amount of DNA in each sample was determined as the ratio of the 2C channel number of the sample and standard times 328 Mbp. The genome size of these I. luminosus males was determined to be 764 ± 7 Mbp (SEM, n = 5). Genome size inference via Kmer spectral analysis of the I. luminosus linked-read data estimated a genome size of 841 Mbp (Appendix 3—figure 2).

HMW DNA (25 µg) 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 µg 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 resulting library was then sequenced on one HiSeqX lane. 408,838,927 paired reads (150 × 150 PE) were produced, corresponding to a genomic coverage of 153x. To evaluate the effect of different Ilumina instruments on data and assembly quality, the library was also sequenced on one HiSeq2500 lane, where 145,250,480 reads (150 × 150 PE) were produced, corresponding to a genomic coverage of 54x. A summary of the library statistics for the genomic sequencing is available in Appendix 4—table 1. The draft genome of I. luminosus (Ilumi1.0) was assembled from the obtained HiSeqX genomic sequencing reads using the Supernova assembler (v1.1.1) (Weisenfeld et al., 2017), on a 40 core 1 TB RAM server at the Whitehead Institute for Biomedical Research. The reported mean molecule size was 12.23 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 long-read based scaffolding was then applied to produce a final assembly Ilumi1.2 (Appendix 3.5.4).

Appendix 3—figure 2

Download asset Open asset

3.5.2 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:97% D:1.3% to C:76.0% D:1.1%. This loss of the endopterygota BUSCOs led us to conclude that the Platyhelminthes annotated scaffolds were authentic scaffolds of I. luminosus, but sequences of Hymenolepis sp. may have been transferred into the I. luminosus genome via horizontal-gene-transfer (HGT). Although Hymenolepis diminuta infects mammals, it also spends a period of its life cycle in intermediate insect hosts, including beetles, as cysticercoids (Center for Disease Control and Prevention, 2017; Sheiman et al., 2006). For a beetle like I. luminosus, which has a extended predatory larval stage, the accidental ingestion and harboring of a Hymenolepis sp. is plausible, potentially enabling HGT between Hymenolepis sp. and I. luminosus over evolutionary timescales.

Appendix 3—figure 3

Download asset Open asset

3.5.4 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). The resulting library was sequenced for 48 hr on a MinION sequencer using a R9.4 flow cell (Part #:FLO-MIN106). Raw trace data was basecalled live within the MinKNOW software (v18.01.6). 824,248 reads (2.4 Gbp; ~1–2x of the I. luminosus genome) were obtained. Reads were mapped to Ilumi1.1 with minimap2 (v2.8-r686-dirty) (Li, 2018) using parameters (-ax map-ont). Inspection of mapped reads with Integrative Genomics Viewer(v2.4.8) (Thorvaldsdóttir et al., 2013) revealed a 17.6 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 non-palindromic 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.

Appendix 3—figure 4

Download asset Open asset

Appendix 3—table 1

Repeat name	Repeat unit length	Repeat unit sequence
Ilumi.complex.repeat.1	~100 bp	TGGTACGAACTATACACGTATACTCAAATCTAATT GTGATACAGCAAAGTAATAATGCAGCATTGTTTGCC GCTCTATACTGCGATTTTATAGTGGT

Appendix 3—figure 5

Download asset Open asset

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

We annotated the coding gene structure of I. luminosus by integrating direct coding gene models produced from the de novo transcriptome (Appendix 3.7.1) and reference guided transcriptome (Appendix 3.7.2), with a lower weighted contribution of ab initio gene predictions, using the Evidence Modeler (EVM) algorithm (v1.1.1) (Haas et al., 2008). First, Augustus (v3.2.2) (Stanke et al., 2006) was trained against Ilumi1.0 with BUSCO (parameters: -l endopterygota_odb9

--long --species tribolium2012). Augustus predictions of Ilumi1.0 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), Photinus pyralis (PPYR_OGS1.0; this report), Aquatica lateralis (AlatOGS1.0; this report), the I. luminosus de novo transcriptome translated peptides, and MAKER blastn/exonerate transcript alignments of the I. luminosus de novo transcriptome transcripts.

We then integrated the ab initio predictions with our de novo and reference guided direct coding gene models, using EVM. In the final version, eight sources of evidence were used for EVM: de novo transcriptome direct coding gene models (Ilumi1.1_Trinity_unstranded-DCGM; weight = 8), reference guided transcriptome direct coding gene models (Ilumi1.1_Stringtie_unstranded-DCGM; weight = 4), MAKER/Augustus ab initio predictions (Ilumi1.1_maker_augustus_ab-initio; weight = 1), protein alignments (P. pyralis, A. lateralis, D. melanogaster, T. castaneum, I. luminosus; weight = 1 each). A custom script (Fallon, 2018a) was used to convert the input MAKER GFF to an EVM compatible GFF format.

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.

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

The mitochondrial genome sequence of I. luminosus was assembled by a targeted sub-assembly 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-sensitive-local ) (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) (Shen et al., 2016) to place the FASTA record break in the AT-rich region, and was submitted to the MITOSv2 mitochondrial genome annotation web server. Small mis-annotations (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).

Appendix 3—figure 6

Download asset Open asset

The level of non-eukaryote contamination of the raw read data for each P. pyralis library was assessed using kraken v1.0 (Wood and Salzberg, 2014) using a dust-masked minikraken database to eliminate comparison with repetitive sequences. Overall contamination levels were low (Appendix 4—table 1), in agreement with a low level of contamination in our final assembly (Appendix 1—figure 9, Appendix 2—figure 2, Appendix 3—figure 3). On average, contamination was 3.5% in the PacBio reads (whole body) and 1.6% in the Illumina reads (only thorax) (Appendix 4—table 1). There was no support for Wolbachia in any of the P. pyralis libraries, with the exception of a single read from a single library which had a kraken hit to Wolbachia. QUAST version 4.3 (Gurevich et al., 2013) was used to calculate genome quality statistics for comparison and optimization of assembly methods (Appendix 4—table 2). BUSCO (v3.0.2) (Simão et al., 2015) was used to estimate the percentage of expected single copy conserved orthologs captured in our assemblies and a subset of previously published beetle genome assemblies (Appendix 4—table 3). The endopterygota_odb9 (metamorphosing insects) BUSCO set was used. The bacteria_odb9 gene set was used to identify potential contaminants by screening contigs and scaffolds for conserved bacterial genes. For genome predictions from beetles, the parameter ‘--species tribolium2012’ was used to improve the BUSCO internal Augustus gene predictions. For Drosophila melanogaster BUSCO genome predictions (Appendix 4—table 3) ‘--species=fly’ was used.

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, octopamine binding, and secreted from the lantern (presumably into the hemolymph of the light organ), could be triggered to release oxygen upon octopamine binding, thereby providing a triggerable O₂ store within the light organ under control of neurotransmitter involved in flash control. As O₂ is believed to be limiting in the adult light reaction, such a release of O₂ could enhance flash intensity or accelerate flash kinetics. Further research is required to test this hypothesis.

Appendix 4—figure 8

Download asset Open asset

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 non-visual 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⁻²⁰) 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.

Appendix 4—figure 9

Download asset Open asset

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² 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-mm-glass 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 µL was obtained. The extracted hemolymph was diluted with 50 µL 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 µL), was diluted with 50 µL MeOH, and air dried. The P. pyralis extracted hemolymph was filtered through a 0.2 μm PFTE filter (Filter Vial, P/No. 15530–100, Thomson Instrument Company), whereas the A. lateralis hemolymph residue was redissolved in 100 µL 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 µL 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 μm 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₂, lyophilized, and the whole body (dry weight: 29.1 mg) was placed in 200 µL 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 μm 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 μm PFTE filter (Filter Vial, P/No. 15530–100, Thomson Instrument Company).

Injections of these filtered extracts (P. pyralis adult male hemolymph 10 µL; A. lateralis adult male hemolymph 5 µL; P. pyralis partial larval body extract 5 µL; A. lateralis whole larval body 5 µL; I. luminosus thorax extract 20 µL) were separated and analyzed using an UltiMate 3000 liquid chromatography system (Thermo Scientific) equipped with a 150 mm C18 Column (Kinetex 2.6 μm silica core shell C18 100 Å pore, P/No. 00F-4462-Y0, Phenomenex, USA) coupled to a Q-Exactive mass spectrometer (Thermo Scientific, USA). 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₂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₂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¹ scan from m/z 120–1250 followed by 1–3 data-dependent MS² scans using HCD fragmentation with a stepped collision energy of 10, 15, 25 normalized collision energy (NCE). Positive mode and negative mode MS¹ and MS² 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 (μA), 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).

Within MZmine2, data were from all five samples on positive mode, and were first cropped to 20 min in order to compare data which was obtained with the same LC gradient parameters. Profile MS¹ data was then converted to centroid mode with the Mass detection module(Parameters: Mass Detector = Exact mass, Noise level = 1.0E4), whereas MS² data was converted to centroid mode with (Noise level = 1.0E1). Ions were built into chromatograms using the Chromatogram Builder module with parameters (min_time_span = 0.10,min_height = 1.0E4, m/z tolerance = 0.001 m/z or five ppm. Chromatograms were then deconvolved using the Chromatogram deconvolution module with parameters (Algorithm = Local Minimum Search, Chromatographic threshold = 5.0%, Search Minimum in RT range = 0.10 min, Minimum relative height = 1%, Minimum absolute height = 1.0E0, Min ratio of peak top/edge = 2, Peak duration range = 0.00–10.00). Isotopic peaks were annotated to their parent features with the Isotopic peaks grouper module with parameters (m/z tolerance = 0.001 or five ppm, Retention time tolerance = 0.2 min, Monotonic shape = yes, Maximum charge = 2, Representative isotope = Most intense). The five peaklists (P. pyralis hemolymph, P. pyralis larval partial body, A. lateralis adult hemolymph, A. lateralis larval whole body, I. luminosus thorax) were then joined and retention time aligned using the RANSAC algorithm with parameters (m/z tolerance = 0.001 or 10 ppm, RT tolerance = 1.0 min, RT tolerance after correction = 0.1 min, RANSAC iterations = 100, Minimum number of points = 5%, Threshold value = 0.5). These aligned peaklists were then gap-filled. Systematic mass accuracy error was determined with the endogenous tryptophan [M + H]⁺ ion (m/z = 205.09, RT = 3.5–4.5 mins), and was measured to be +0.6 ppm,+9.9 ppm,+1.6 ppm,+1.1 ppm, and +0.6 ppm, for P. pyralis adult hemolymph, P. pyralis partial larval body extract, A. lateralis adult hemolymph, A. lateralis larval body extract, and I. luminosus thorax extract, respectively.

Appendix 4—figure 10

Download asset Open asset

Appendix 4—figure 11

Download asset Open asset

Appendix 4—figure 12

Download asset Open asset

Appendix 4—figure 13

Download asset Open asset

Appendix 4—figure 14

Download asset Open asset

Appendix 4—figure 15

Download asset Open asset

Appendix 4—figure 16

Download asset Open asset

Appendix 4—figure 17

Download asset Open asset

Appendix 4—figure 18

Download asset Open asset

Appendix 4—figure 19

Download asset Open asset

4.6.5 MS² similarity search for P. pyralis lucibufagins

We first performed a MS² similarity search within P. pyralis adult hemolymph for ions that showed a similar MS² spectra to the MS² spectra arising from the diacetylated lucibufagin [M + H]⁺ ion from the same run ([M + H]⁺ m/z 533.2385, RT = 15.10 mins) (Appendix 4—figure 20). This search was performed through the MS² similarity search module of MZmine2 (v2.30) with parameters (m/z tolerance: 0.0004 m/z or 1 PPM; minimum # of ions to report: 3). This MS² similarity search revealed nine putative lucibufagin isomers with highly similar MS² spectra (Appendix 4—figure 21), which expanded to 17 putative lucibufagin isomers when considering features without MS² spectra, but with identical exact masses and close retention times (ΔRT <2 min) to the previously identified 9 (Appendix 4—table 7). Chemical formula prediction was assigned to each precursor ion using the Chemical formula search module of MZmine2, whereas chemical formula predictions for product ions was performed within MZmine2 using SIRUIS (v3.5.1) (Böcker et al., 2009). The structural identity of the nine putative lucibufagins detected via the MS² spectra similarity search was easily interpreted in light that the different chemical formula represented the core lucibufagins that had undergone acetylation (COCH₃) or propylation (COCH₂CH₃), in different combinations. Notably, the most substituted isomers, dipropylated lucibufagin ([M + H]⁺ m/z 561.2695, RT = 19.54 mins) were close to the edge of the cropped data (20 min), thus it may be possible that more highly substituted lucibufagins with a longer retention times are present, but not detected in the current analysis.

We then performed a MS² similarity search within P. pyralis partial body extract for ions that showed a MS² spectra similar to that of the dipropylated lucibufagin [M + H]⁺ ion from the same run ([M + H]⁺ m/z 561.2738, RT = 19.53). This search was performed through the MS² similarity search module of MZmine2 (v2.30) with parameters (m/z tolerance: 0.0004 m/z or 1 PPM; minimum # of ions to report: 5). This MS² similarity search revealed 14 putative lucibufagin isomers with highly similar MS² spectra (Appendix 4—table 7). Complexes and fragments were manually removed from the analysis. Comparison of the theoretical and observed exact mass indicated that this experimental run had an unusual degree of systematic m/z error, of ~+10 ppm. After manual correction m/z, chemical formula prediction revealed a several putative lucibufagins of unknown structure with nitrogen in their chemical formula, suggesting that the nitrogen containing lucibufagins reported by by Gronquist and colleagues from Lucidota atra (Gronquist et al., 2005) may be present in P. pyralis larvae.

Appendix 4—figure 20

Download asset Open asset

Appendix 4—figure 21

Download asset Open asset

Appendix 4—table 7

Assigned ion identity	Ion type	Chemical formula	Expected M/z	Measured M/z	M/z error* (ppm)	Retention time (mins)	Feature area (arb)
Core lucibufagin isomer 1	[M + H]+	C24H33O8	449.2175	449.2171	−0.89	7.9	6.7E + 05
Core lucibufagin isomer 2	""	“”	“”	""	“”	9.3	1.1E + 07
Monoacetylated lucibufagin isomer 1	""	C26H35O9	491.2281	491.2277	−0.81	10.2	4.2E + 07
Core lucibufagin isomer 3	""	C24H33O8	449.2175	449.2171	−0.89	10.8	1.7E + 07
Monoacetylated lucibufagin isomer 2	""	C26H35O9	491.2281	491.2277	−0.81	11.4	1.1E + 06
Monoacetylated lucibufagin isomer 3	""	“”	“”	""	“”	11.9	1.8E + 07
Monoacetylated lucibufagin isomer 4	""	“”	“”	""	“”	13.0	2.7E + 08
Monoacetylated lucibufagin isomer 5	""	“”	“”	""	“”	13.2	6.0E + 07
Monoacetylated lucibufagin isomer 6	""	“”	“”	""	“”	14.5	6.2E + 06
Diacetylated lucibufagin isomer 1	""	C28H37O10	533.2387	533.2385	−0.37	15.1	4.0E + 09
Diacetylated lucibufagin isomer 2	""	“”	“”	""	“”	15.4	1.9E + 09
Monoacetylated, mono propylated lucibufagin isomer 1	""	C29H39O10	547.2543	547.2542	−0.18	17.0	1.5E + 07
Monoacetylated, mono propylated lucibufagin isomer 2	""	“”	“”	""	“”	17.4	2.8E + 08
Monoacetylated, mono propylated lucibufagin isomer 3	""	“”	“”	""	“”	17.7	1.2E + 08
Dipropylated lucibufagin isomer 1	""	C30H41O10	561.2700	561.2695	−0.89	18.9	1.4E + 08
Dipropylated lucibufagin isomer 2	""	“”	“”	""	“”	19.5	3.9E + 07
Dipropylated lucibufagin isomer 3	""	“”	“”	""	“”	19.8	1.8E + 08

Appendix 4—table 8

Assigned ion identity	Ion type	Chemical formula	Expected m/z	Measured m/z	m/z error (ppm)	Retention time (mins)	Feature area (arb)
Core lucibufagin isomer 2	[M + H]+	C24H33O8	449.2175	449.2215	+8.9	9.15	8.5E + 06
Monoacetylated lucibufagin isomer 1	“”	C26H35O9	491.2277	491.2326	+9.9	10.04	1.2E + 07
Unknown	unknown	C28H39O10*	535.2543*	535.2592	+9.1*	12.40	1.6E + 07
Unknown	unknown	C24H38NO6*	436.2695*	436.2735	+9.1*	13.30	2.2E + 07
Unknown	unknown	C27H45N2O8*	525.3173*	525.3221	+9.1*	13.35	1.3E + 08
Unknown	unknown	C24H40NO7*	454.2799*	454.2840	+9.1*	13.73	1.3E + 07
Diacetylated lucibufagin isomer 1	[M + H]+	C28H37O10	533.2387	533.2426	+7.3	14.93	1.7E + 09
Diacetylated lucibufagin isomer 2	[M + H]+	“”	“”	533.2426	+7.3	15.16	3.5E + 08
Unknown	Unknown	C29H46NO8*	536.3216*	536.3256	+7.3*	16.57	4.1E + 07
Unknown	Unknown	Unknown**	563.2854*	563.2896	+7.3*	16.80	1.3E + 07
Unknown	Unknown	C26H31O7	455.2056	455.2097	+9.1*	17.22	5.8E + 07
Dipropylated lucibufagin isomer 3	Unknown	C30H41O10	561.2700	561.2738	+6.7	19.53	2.0E + 09
Dipropylated lucibufagin isomer 4	Unknown	C30H41O10	561.2700	561.2738	+6.7	19.82	2.2E + 08

Appendix 4—table 9

Chemical formula	Predicted exact mass	Exact mass adjusted to P. pyralis hemolymph data (+0.6 ppm)	Exact mass adjusted to P. pyralis partial larval body data (+9.9 ppm)	Exact mass adjusted to A. lateralis hemolymph data (+1.6 ppm)	Exact mass adjusted to A. lateralis larval body data (+1.1 ppm)	Exact mass adjusted to I. luminosus thorax data (+0.6 ppm)
C24H33O8	449.2175	449.2178	449.2219	449.2182	449.2180	449.2178
C24H38NO6*	436.2699	436.2702	436.2742	436.2706	436.2704	436.2702
C24H40NO7*	454.2804	454.2807	454.2849	454.2811	454.2809	454.2807
C26H31O7	455.2069	455.2072	455.2114	455.2076	455.2074	455.2072
C26H35O9	491.2281	491.2284	491.2330	491.2289	491.2286	491.2284
C27H45N2O8*	525.3175	525.3178	525.3227	525.3183	525.3181	525.3178
C28H37O10	533.2386	533.2389	533.2439	533.2395	533.2392	533.2389
C28H39O10*	535.2543	535.2546	535.2596	535.2552	535.2549	535.2546
C29H39O10	547.2543	547.2546	547.2597	547.2552	547.2549	547.2546
C29H46NO8*	536.3223	536.3226	536.3276	536.3232	536.3229	536.3226
C30H41O10	561.2699	561.2702	561.2755	561.2708	561.2705	561.2702

1. *=Chemical formula assigned for structurally unclear putative lucibufagins

The complete genome of the molicute (Phylum: Tenericutes) Entomoplasma luminosum subsp. pyralis was constructed by a long-read metagenomic sequencing and assembly approach from the P. pyralis PacBio data. First, BUSCO v.3 with the bacterial BUSCO set was used to identify those contigs from the PacBio only Canu assembly (Ppyr0.1-PB) which contained conserved bacterial genes. A single 1.04 Mbp contig with 73 bacterial BUSCO genes was the only contig identified with more than 1 BUSCO hit. Inspection of the Canu produced assembly graph with Bandage v0.8.1 (Wick et al., 2015), revealed that the contig had a circular assembly path. BLASTN alignment of the contig to the NCBI nt database indicated that this contig had a high degree of similarity to annotated Mycoplasmal genomes. Together this data suggested that this contig represented a complete Mycoplasmal genome. Polishing of the contig was performed by mapping and PacBio consensus-calling using SMRTPortal v2.3.0.140893 with the ‘RS_Resequencing.1’ protocol with default parameters. The median coverage was ~50x. The resulting consensus sequence was restarted with seqkit (Shen et al., 2016) to place the FASTA record junction 180˚ across the circular chromosome, and reentered into the polishing process to enable efficient mapping across the circular junction. This mapping, consensus calling, and rotation process was repeated three times total, after which no additional nucleotide changes occurred. The genome was ‘restarted’ with seqkit such that the FASTA start position began between the ribosomal RNAs, and annotation was conducted through NCBI using their prokaryotic gene annotation pipeline (PGAP). Analysis with BUSCO v.3 of the peptides produced from the aforementioned genome annotation indicated that 89.8% of expected Tenericutes single-copy conserved orthologs were captured in the annotation (C:89.8%[S:89.8%,D:0.0%], F:2.4%, M:7.8%, n:166). Comparison of the predicted 16S RNA gene sequence to the NCBI 16S RNA gene database indicated that this gene had 99% identity to the E. luminosum 16S sequence (ATCC 49195 - formerly Mycoplasma luminosum; NCBI Assembly ID ASM52685v1) (Kyrpides et al., 2014; Williamson et al., 1990), leading to our description of this genome as the genome of Entomoplasma luminosum subspecies (subsp.) pyralis. Protein overlap comparisons using the OrthoFinder pipeline (v1.1.10) (Emms and Kelly, 2015) between our predicted protein geneset for E. luminosum var. pyralis and the protein geneset of Entomoplasma luminosum (ATCC 49195 - formerly M. luminosum; NCBI Assembly ID ASM52685v1), indicated that 94% (670/709) of the previously annotated E. luminosum proteins are present in our genome of E. luminosum subsp. pyralis.

The complete mitochondrial genome of the dipteran parasatoid Apocephalus antennatus, first detected via BLASTN of mtDNAs as a concatemerized sequence in the Canu PacBio only assembly (Ppyr0.1-PB) was constructed in full by a long-read metagenomic sequencing and assembly approach. First, PacBio reads were mapped to the NCBI set of mitochondrial genomes concatenated with the P. pyralis mitochondrial genome assembly reported in this manuscript (NCBI accession KY778696.1), using GraphMap v0.5.2 with parameters ‘align -C -t 4 -P’. Of the mitochondrially mapped reads (45949 reads), 98% (45267 reads) were partitioned to the P. pyralis mtDNA. The next most abundant category at 1.1% (531 reads), was partitioned to the mtDNA of the Phorid fly Megaselia scalaris (NCBI accession: KF974742.1). The next most abundant category at 0.11% (53 reads) was partitioned to the mitochondrion of the Red algae Galdieria sulphuraria (NCBI accession: NC_024666.1). The reads were then split into three partitions: P. pyralis mapping, M. scalaris mapping, and other, and input into Canu (v1.6+44) (Koren et al., 2017) for assembly. Each partitioned assembly by Canu produced a single circular contig, notably the ‘other’ and Megaselia partitions produced highly similar sequences, whereas the P. pyralis partition produced a circular sequence that was highly similar to the P. pyralis mtDNA. We inspected the M. scalaris partition further as it was produced with more reads. Notably, although an inspection of the contig was circular, and showed a high degree of similarity upon blastn to the M. scalaris mtDNA, the contig was ~2x larger than expected (29,821 bp). An analysis of contig’s self-complementarity with Gepard (v1.40) (Krumsiek et al., 2007), indicated that this contig had 2x tandem repetitive regions, and was duplicated overall twice. Similarly, the. GFA output of Canu noted an overlap of 29,821, indicating that the assembler was unable to determine an appropriate overlap, other than the entire contig. Manual trimming of the contig to the correct size, 180˚ restarting with seqkit, and polishing using SMRTPortal v2.3.0.140893 with the ‘RS_Resequencing.1’ protocol with default parameters, followed by 180˚ seqkit ‘restarting’, followed by another round of polishing, produced the final mtDNA (18,674 bp; Appendix 5—figure 1). This mtDNA was taxonomically identified in a separate analysis to originate from A. antennatus (Appendix 5.3). 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).

Appendix 5—figure 1

Download asset Open asset

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.

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). They are multipartite linear ssRNA negative strand viruses, conformed by five genome segments generating a ca. 10.8 Kbp total RNA genome. Genome segments one through three (ca. 2.3–2.5 Kbp long) encode a heterotrimeric viral polymerase constituted by subunit Polymerase Basic protein 1 - PB1 (PpyrOMLV1: 801 aa, 91 kDA; PpyrOMLV2: 802 aa, 91.2 kDA), Polymerase Basic protein 2 - PB2 (PpyrOMLV1: 804 aa, 92.6 kDA; PpyrOMLV2: 801 aa, 92.4 kDA) and Polymerase Acid protein - PA (PpyrOMLV1: 754 aa, 86.6 kDA; PpyrOMLV2: 762 aa, 87.9 kDA). PpyrOMLV1 and PpyrOMLV2 PB1 present a Flu_PB1 functional domain (Pfam: pfam00602; PpyrOMLV1: interval = 49–741, e-value = 2.93e-69; PpyrOMLV2: interval = 49–763, e-value = 1.42e-62) which is the RNA-directed RNA polymerase catalytic subunit, responsible for replication and transcription of virus RNA segments, with two nucleotide-binding GTP domains. PpyrOMLV1 and PpyrOMLV2 PB2 present a typical Flu_PB2 functional domain (Pfam: pfam00604; PpyrOMLV1: interval = 26–421, e-value = 5.10e-13; PpyrOMLV2: interval = 1–692, e-value = 1.57e-11) which is involved in 5' end cap RNA structure recognition and binding to further initiate virus transcription. PpyrOMLV1 and PpyrOMLV2 PA subunits share a characteristic Flu_PA domain (Pfam: pfam00603; PpyrOMLV1: interval = 122–727, e-value = 3.73e-07; PpyrOMLV2: interval = 117–732, e-value = 5.63e-10) involved in viral endonuclease activity, necessary for the cap-snatching process (Guilligay et al., 2014). Genome segment four (1.6 Kbp size) encodes a Hemaglutinin protein – HA (PpyrOMLV1: 526 aa, 59.7 kDA; PpyrOMLV2: 525 aa, 58.6 kDA) presenting a Baculo_gp64 domain (Pfam: pfam03273; PpyrOMLV1: interval = 108–462, e-value = 2.16e-15; PpyrOMLV2: interval = 42–460, e-value = 1.66e-23), associated with the gp64 glycoprotein from baculovirus as well as other viruses, such as Thogotovirus (Orthomyxoviridae - OMV) which was postulated to be related to the arthropod-borne nature of these specific Orthomyxoviruses. In addition, HA as expected, presents an N-terminal signal domain, a C terminal transmembrane domain, and a putative glycosylation site. Lastly, genome segment five (ca. 1.8 Kbp size) encodes a putative nucleocapsid protein – NP (PpyrOMLV1: 562 aa, 62.3 kDA; PpyrOMLV2: 528 aa, 58.5 kDA) with a Flu_NP structural domain (Pfam: pfam00506; PpyrOMLV1: interval = 145–322, e-value = 1.32e-01; PpyrOMLV2: interval = 94–459, e-value = 1.47e-04) this single-strand RNA-binding protein is associated to encapsidation of the virus genome for the purposes of RNA transcription, replication and packaging (Appendix 5—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′ 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 (Reich et al., 2017). 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, incorporating virus glycoproteins and eventually non-glycosylated proteins (one or two in number). Typical virion surface glycoprotein projections are 10–14 nm in length and 4–6 nm in diameter. The virus genome is multisegmented, has a helical-like symmetry, consisting of different size ribonucleoproteins (RNP), 50–150 nm in length. Influenza RNPs can perform either replication or transcription of the same template. Virions of each genus contain different numbers of linear ssRNA (-) genome segments (King et al., 2011). Influenza A virus (FLUAV), influenza B virus (FLUBV) and infectious salmon anemia virus (ISAV) are conformed of eight segments. Influenza C virus (FLUCV), Influenza D virus (FLUDV) and Dhori virus (DHOV) have seven segments. Thogoto virus (THOV) and Quaranfil virus (QUAV) have six segments. Johnston Atoll virus (JAV) genome is still incomplete, and only two segments have been described. Segment lengths range from 736 to 2396 nt. Genome size ranges from 10.0 to 14.6 Kbp (King et al., 2011). As described previously, every OMV RNA segment possess conserved and partially complementary 5′- and 3′-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 non-glycosylated 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 sialyl-galactosyl 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 (Reich et al., 2017; 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′ end. Virus mRNAs are polyadenylated at the 3′ 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 re-evaluate 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).

Appendix 5—figure 2

Download asset Open asset

Appendix 5—figure 3

Download asset Open asset

To gain insights on the potential shared evolutionary history of P. pyralis and the IOMV PpyrOMLV1 and 2, we examined our assembly of P. pyralis for putative signatures or paleovirological traces (Ballinger et al., 2014; Metegnier et al., 2015; Feschotte and Gilbert, 2012) that would indicate ancestral integration of virus related sequences into the firefly host. Remarkably, we found Endogenous virus-like Elements (EVEs) (Katzourakis and Gifford, 2010), sharing significant sequence identity with most PpyrOMLV1 and 2 genome segments, spread along four P. pyralis linkage-groups. Virus integration into host genomes is a frequent event derived from reverse transcribing RNA viruses (Retroviridae). Retroviruses are the only animal viruses that depend on integration into the genome of the host cell as an obligate step in their replication strategy (Temin, 1985). Viral infection of germ line cells may lead to viral gene fragments or genomes becoming integrated into host chromosomes and subsequently inherited as host genes.

Animal genomes are paved by retrovirus insertions (Bushman et al., 2005). These insertions, which are eventually eliminated from the host gene pool within a few generations, and may, in some cases, increase in frequency, and ultimately reach fixation. This fixation in the host species can be mediated by drift or positive selection, depending on their selective value. On the other hand, genomic integration of non-retroviral viruses, such as PpyrOMLV1 and 2, is less common. Viruses with a life cycle characterized by no DNA stage, such as OMV, do not encode a reverse transcriptase or integrase, thus are not retro transcribed nor integrated into the host genome. However, exceptionally and recently, several non-retroviral sequences have been identified on animal genomes; these insertions have been usually associated with the transposable elements machinery of the host, which provided a means to genome integration (Gilbert and Cordaux, 2017; Palatini et al., 2017). Interestingly, when we screened our P. pyralis genome assembly Ppyr1.2 by BLASTX searches (E-value <1e10⁻⁶) of PpyrOMLV1 and 2 genome segments, we identified several genome regions that could be defined as Firefly EVEs, which we termed FEVEs (Appendix 5—figure 2J; Appendix 5—table 5-8). We found 30 OMV related FEVEs, which were mostly found in linkage group one (LG1, 83% of pinpointed FEVEs). The majority of the detected FEVEs shared sequence identity to the PB1 encoding region of genome segment one of PpyrOMLV1 and 2 (ca. 46% of FEVEs; Appendix 5—table 5), followed by NP encoding genome segment five (ca. 33% of detected FEVEs; Appendix 5—table 8). In addition we identified four FEVEs related to genome segment three (PA region; Appendix 5—table 7) and two FEVEs associated to genome segment two (PB2 encoding region; Appendix 5—table 6). We found no evidence of FEVEs related to the hemagglutinin coding genome segment four (HA) via BLASTX. The detected P. pyralis FEVEs represented truncated fragments of virus like sequences, generally presenting frameshift mutations, early termination codons, lacking start codons, and sharing diverse mutations that altered the potential translation of eventual gene products. FEVEs shared sequence similarity to the coding sequence of specific genome segments of the cognate FOLMV. We generated best/longest translation products of the corresponding FEVEs, which presented an average length of ca. 21.86% of the corresponding PpyrOMLV genome segment encoding gene region (Appendix 5—table 5-5.5.5), and an average pairwise identity to the FOLMV virus protein of 55.08%. Nevertheless, we were able to 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 EVE-dependent 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.