Genomic and transcriptomic resources for assassin flies including the complete genome sequence of Proctacanthus coquilletti (Insecta: Diptera: Asilidae) and 16 representative transcriptomes

Genomes available for Diptera

To date, out of the 160,000 species of Diptera (Pape, Blagoderov & Mostovski, 2011), complete genomes have been sequenced for 100 species (NCBI as of 04 October 2016). These represent 47% of the insect genomes available at NCBI and are concentrated within the earliest radiation of Diptera including the medically important mosquitoes (Culicidae) and sand flies (Psychodidae) and the higher flies including the model organism Drosophila and medically important Glossina tsetse flies (Fig. 1).

A recent study by Vicoso & Bachtrog (2015) added some 37 low-coverage genomes for a study on the sex chromosomes of Diptera. While the genome sequencing in this publication was not intended to add draft genomes, the genomes are spread across the Diptera Tree of Life (Fig. 1) and Vicoso & Bachtrog (2015) added two low-coverage (approximately 12×) genomes for Orthorrhapha in the center of fly evolution, i.e., the soldier fly Hermetia illucens (Stratiomyomorpha, estimated genome size = 1.3 Gbp, N50 = 2,778 bp) and the assassin fly Holcocephala fusca (Asiloidea, estimated genome size = 673 Mbp, N50 = 4,591 bp).

Aedes, Anopheles, and Culex mosquitoes and Drosophila vinegar flies shared a common ancestor some 240 Million years ago (mya) (Wiegmann et al., 2011). The most recent common ancestor of mosquitoes and Orthorrhapha likewise lived 240 mya and that of Drosophila and Orthorrhapha some 200 mya. Filling a gap in the center of the Diptera Tree of Life (Fig. 1) by providing data on novel, high-quality draft genomes from within Orthorrhapha and Asiloidea will open the opportunity to more meaningfully compare genomes across Diptera. Furthermore, the genomic resources provided here will advance the study of evolutionary history, life history, and the search for the venom of assassin flies.

Specimen source

Adult flies were hand-netted either directly from their resting/perching sites (Apioceridae, Asilidae, Bombyliidae, and Mydidae) or from within a Malaise Trap (Tabanidae) and kept alive in individual vials. They were identified to species, assigned unique identifiers, and either preserved in RNAlater (specimen cut open and placed directly in RNAlater) or liquid N₂ (specimen alive in individual vial dropped in dry shipper containing liquid N₂). RNAlater vials were emptied of any liquid before being placed in liquid N₂-filled tanks in the NMNH Biorepository where all specimens are stored and accessible by their unique specimen identifier (Table 1).

RNA-Seq

Total RNA was extracted from specimens preserved in RNAlater or in liquid N₂ (see Table 1). A single specimen was used for each extraction. Muscular tissue was extracted from the thorax and cryogenically ground using CryoMill (Retsch, Haan, Germany). Total RNA was isolated using the TRI Reagent Protocol (Sigma-Aldrich, St. Louis, MO, USA) with overnight precipitation, and then quantified using Epoch Microplate Spectrophotometer and Gen5 software (both BioTek, Winooski, VT, USA). For the specimens sequenced using Ion Torrent, the isolation of mRNA was carried out using DynaBeads mRNA DIRECT Kit, and Ion Total RNA-Seq Kit (v2) for Whole Transcriptome Libraries (Thermo Fisher Scientific) was used for library preparation. The BluePippin System (Sage Science, Beverly, MA, USA) was used for selecting fragments of 170–350 bp. For the specimens sequenced using Illumina MiSeq and HiSeq2000, the isolation of mRNA and construction of stranded mRNA-Seq libraries were carried out using KAPA Stranded mRNA-Seq Kit (Kapa Biosystems, Boston, MA, USA) and NEBNext Multiplex Oligos (New England BioLabs, Ipswich, MA, USA). Library fragment size distribution was assessed using High Sensitivity D1000 ScreenTape System (Agilent, Waldbronn, Germany), and the BluePippin System was used for selecting fragments of 180–440 bp. After size selection, a sample of each library was quantified using the KAPA Library Quantification Kit for Illumina platforms, and pooled to 5 nM total concentration for sequencing.

RNA-Seq bioinformatics workflow is shown in Fig. 2. Raw data as well as assembled transcripts were screened for contamination with KRAKEN (Wood & Salzberg, 2014). RNA-Seq reads were trimmed with Trimmomatic (Bolger, Lohse & Usadel, 2014) and transcripts were assembled in Trinity (Grabherr et al., 2011). Transcriptome “completeness” was estimated using BUSCO (v2.0beta, Simão et al. (2015)) with the “Endopterygota” lineage specific set of 2,442 loci and the -m tran setting. BUSCO assesses completeness with near-universal single copy orthologs selected from OrthoDB (Kriventseva et al., 2015). The 16 sets of assembled transcripts were translated with Transdecoder (Haas & Papanicolao, 2015) under default parameters. Peptides were filtered for redundancy using CD-Hit (v4.6.1, Fu et al. (2012)) specifying a 95% similarity threshold. The Trinotate workflow was used for transcript annotation (Grabherr et al., 2011). Trinotate uses evidence from BLASTx, BLASTp (Altschul et al., 1990), PFAM (Punta et al., 2012), and HMMER (Finn, Clements & Eddy, 2011) to assign GO terms (Ashburner et al., 2000) to transcripts. Venn diagrams showing overlapping sets of GO-terms were generated at: http://bioinformatics.psb.ugent.be/webtools/Venn/ (see Fig. 3).

Genome sequencing

Genomic DNA was extracted from thoracic muscular tissue and legs of a single specimen of Proctacanthus coquilletti preserved in liquid N₂, using the DNEasy DNA Extraction kit (Qiagen, Hilden, Germany). The sample was quantified using Epoch Microplate Spectrophotometer and Gen5 software, and subsequently pooled to 50 ng/µL concentration. Sequencing took place at Johns Hopkins University Genetic Resources Core Facility’s High Throughput Sequencing Center. A PCR-free library was generated and two lanes of Illumina HiSeq2500 were sequenced to satisfy DISCOVAR recommendations.

Genome sequencing bioinformatics workflow is shown in Fig. 2. Genome size, heterozygosity, and repeat content were estimated with raw reads using GenomeScope (Sedlazeck, Nattestad & Schatz, 2016), which uses a kmer histogram generated by JELLYFISH (Marçais & Kingsford, 2011). Raw data as well as assembled contigs were screened for contamination with KRAKEN. Blobtools/Blobology (Kumar et al., 2013) was also used to assess contamination. Sequences were assembled using DISCOVAR de novo (Jaffe, 2015) and w2rap-contigger (Clavijo, 2016) with a kmer size of 200 and 260. w2rap-contigger provides performance improvements on DISCOVAR de novo, which is no longer being actively developed. Some scaffolding was performed as with 2 × 250 bp reads there is some space between overlap and DISCOVAR de novo and w2rap-contigger perform scaffolding internally as shown in Fig. 4. Genome completeness was estimated using BUSCO with the “Endopterygota” lineage specific set of loci and the -m genome setting.

Gene prediction was performed using MAKER (Cantarel et al., 2008), which uses RepeatMasker (Smit, Hubley & Green, 2013), Augustus (Stanke & Waack, 2003), BLAST (Altschul et al., 1990), Exonerate (Slater & Birney, 2005), and SNAP (Korf, 2004). Contigs shorter than 2 kbp were not annotated, as they are too short to produce high-quality evidence. The maximum intron size was set as 10 kbp, which is recommended based on Drosophila intron sizes. The Augustus model species used was “fly” and Drosophila melanogaster Repbase libraries (Bao, Kojima & Kohany, 2015) were used for RepeatMasker.

Blobplots

A blobplot was created using blobtools (Kumar et al., 2013, v 0.9.19.5). Prior to generating the blobplot, two steps need to be taken: (1) the raw reads need to be mapped back to the genome to generate an estimate of sequencing coverage and (2) a taxon assignment for each contig needs to be generated by querying the NCBI nucleotide database. Raw reads were mapped using Bowtie 2 (Langmead & Salzberg, 2012, v 2.2.9) and taxon assignments generated using megablast (Altschul et al., 1990; Zhang et al., 2000).

Phylogenetic trees

Orthology detection was conducted in OMA standalone (v1.0.6, http://omabrowser.org/standalone/) using peptides processed in Transdecoder and CD-Hit. Amino acid alignments on individual resulting orthologs were conducted in MAFFT (Katoh, Asimenos & Toh, 2009). Phylogenetic model selection was performed with PartitionFinderProtein (Lanfear et al., 2012). Gene trees and trees of concatenated, partitioned, data were built in RAxML (raxmlHPC-PTHREADS-SSE3, Stamatakis (2014)). Best tree searches were run 100 times each and rapid bootstrapping was run under the AutoMRE option. ASTRAL (Mirarab et al., 2014) was used to generate a species tree based on individual gene tree topologies.

Fossil calibrations

Seven fossils ranging in age from 112–45 Million years old (myo) were used to calibrate the time-tree analysis. The maximum age of the root was set to 195 million years (my), an age proposed for the most recent common ancestor of Tabanidae and Asilidae (Wiegmann et al., 2011). MCMCtree, part of the PAML package (Yang, 2007), was used to generate a time-calibrated phylogeny based on the best-scoring RAxML tree.

Sequence, genome, and analysis data

The raw and assembled sequence data can be accessed under NCBI Umbrella BioProject PRJNA345052. Individual BioProject and BioSample accession numbers are provided in Table 2. The Proctacanthus coquilletti de novo genome assembly (w2rap-contigger 200 kmer) can be accessed under NCBI WGS MNCL00000000 and the Genomescope results at: http://qb.cshl.edu/genomescope/analysis.php?code=TRpKdSHytjlB1vBGsPne. Digital copies of visualizations, alignments, and phylogenetic trees can be accessed under a Figshare Collection (doi: 10.6084/m9.figshare.c.3521787, Table S1).

RNA-Seq

RNA-Seq results are shown in Table 3. It is clear that the Ion Torrent platform was the far inferior sequencing platform in terms of total reads, transcripts, and BUSCO recovered loci. We include these four species in the phylogenomic analyses discussed later in spite of the poor wet-lab results because they still provide sufficient homology to generate a phylogenetic hypothesis even though the number of orthologs is small. The MiSeq and HiSeq results are quite comparable, showing that the larger number of reads generated by the HiSeq is not necessary to achieve strong BUSCO results. Four taxa were pooled on each HiSeq lane, and produced approximately three times the raw data of each MiSeq run. For all samples, data were gathered based on a single life-stage, and therefore do not represent the complete transcriptome, which would require larvae and multiple replications. Larvae are almost impossible to find, making this goal unlikely to be achieved. However, for a number of species, P. albiceps, P. coquilletti, T. atricapillus, T. discus, E. pinguis, in particular, very high-quality sets of transcripts were produced that represent the overwhelming majority of BUSCOs.

The number of GO terms (duplicates removed) in common (overlapping) or unique (not overlapping) for four separate subsets of taxa are summarized in Fig. 3. Figure 3A shows representatives from the major lineages included, i.e., , Apiocera (Apioceridae), Ectyphus (Mydidae), Fidena (Tabanidae), Proctacanthus and Scleropogon (separate Asilidae clades). Figure 3B shows just one representative from each Asiloid family, Fig. 3C the “outgroup” taxa (Tabanidae and Bombyliidae), and Fig. 3D Asilidae exclusive of the Asilinae clade. Considering that the most recent common ancestor of horse flies (Tabanidae) and asiloid flies (Apioceridae, Asilidae, and Mydidae) existed some 195 mya, the GO term overlap among these taxa is quite large with 12,401 terms in common to the best-sequenced taxa (Fig. 3A). Interestingly, the species selected for genome sequencing, Proctacanthus coquilletti, has some 1,149 unique GO terms (Fig. 3A). Among the asiloid flies, some 13,185 shared GO terms were found (Fig. 3B) and here again Proctacanthus coquilletti shows a very high number of unique terms (1,767)—the highest in our comparative study. Similarly high unique GO term numbers are found in the sequenced horse flies (Fig. 3C), which might suggest that blood-sucking and predatory flies have a higher number of unique terms compared to nectar-eating flies such as Apioceridae and Mydidae (Fig. 3B) or Bombyliidae (Fig. 3C). The shared GO terms in Fig. 3D for a clade of Asilidae are much lower, which is most likely caused by the reduced number of Ion Torrent reads for Lasiopogon and Nicocles.

A. parkeri, the only species of Apioceridae sequenced, produced more than two times the number of transcripts than any other taxon (Table 3). A. parkeri also has the most raw reads, but for other species the number of reads does not show the same relationship to the number of transcripts. A. parkeri did not have the highest BUSCO score (i.e., did not produce the most complete single copy orthologs of all the sequenced transcriptomes), but it produced a large number of transcripts with assigned GO terms (14,571), which is the 7th highest number. A. parkeri does not have an unusual number of unique GO assignments (compared to the other taxa we sampled, Figs. 3A–3B). F. pseudoaurimaculata and T. discus, the two species of Tabanidae included, have higher number of transcripts and GO terms. With fewer total transcripts, Mydidae (Ectyphus pinguis, Mydas clavatus, and Messiasia californica, Table 3), the sister group to Apioceridae, reach a high number of GO terms that are with one exception higher than for Apiocera. Interestingly, the only included fly with parasitoid larvae, the bee fly Thevenetimyia, has by far the lowest number of GO terms (Table 3), however, this transcriptome was sequenced on the Ion Torrent platform. As more species from these taxa are sequenced for transcriptomes and genomes we will gain the ability to investigate whether this pattern will hold and why it might be. Horse flies, mydas flies, and apiocerid flies have very different life histories than assassin flies. Almost all females of horse flies are blood-feeders as adults, while males feed on nectar (Lessard et al., 2013; Morita et al., 2015). Adult flies of Apioceridae and Mydidae are nectar- or honeydew-feeders (Norris, 1936; Paramonov, 1953; Cazier, 1982; Wharton, 1982), if they feed at all.

As a quality check, we ran KRAKEN on all assembled transcripts, which for each species resulted in approximately 0.5% of transcripts having any kmer match to a database of all finished RefSeq bacterial, archaeal, and viral genomes at NCBI. This left us with confidence that we are not including contaminants in our numbers in Table 3.

Genome sequencing

One HiSeq2500 flow cell (2 lanes) produced 382,575,358 reads. Pre-assembly assessment of kmer distributions in JELLYFISH to produce a histogram, which is then interpreted by GenomeScope, is summarized in Table 4. GenomeScope provides an estimate of genome size of just under 200 Mbp, a very low rate of heterozygosity (appproximately 0.47%), a small percentage of repeats (approximately 15%), and a very low error rate (0.82%). Before submitting P. coquilletti for sequencing, we did not have a reliable estimate for any of these parameters, and perhaps could have been successful with a single lane of sequencing, since the genome structure does not appear particularly challenging when compared to many other insect genomes. The low-coverage Holcocephala fusca genome (Vicoso & Bachtrog, 2015) is reported to have a genome size of 673 Mbp, more than three times larger than we estimate for P. coquilletti. The contig N50 value for H. fusca is only 1,778 bp, however, making it a less reliable estimate than the one presented here.

Another positive point about assassin flies beyond their relative genomic simplicity is the large thoracic muscle mass that can be used for DNA and RNA extraction. The gut did not have to be included in the extraction to produce enough DNA, which for a PCR-free library is substantial (3 µg). Since the gut is the source of the most obvious contaminants (meals and gut microbiome), this can be an important factor for those sequencing insect genomes.

Genome assembly and annotation

Assembly statistics from DISCOVAR de novo and w2rap-contigger are summarized in Table 5. Our final assembly was produced by the w2rap-contigger 200 kmer assembly. One of the reasons we chose to try w2rap-contigger was that DISCOVAR produced a much larger than expected genome size estimate. We realized that a set of sequences that are appended to the DISCOVAR a.lines.fasta file and are equal to or smaller in length than raw reads and do not represent valid contigs was appended erroneously. DISCOVAR’s own N50 calculator ignores these sequences, and it was only by using our own metadata parser that we found a vastly different N50 value and decided to investigate the assembly output further. We feel that the PCR-free DISCOVAR recipe has produced a very high quality draft genome, particularly given that it is based on a single paired-end library. A plot summarizing the contigs and BUSCO results utilizing the www.lepbase.org assembly statistics tools (Challis, 2016) is shown in Fig. 4 and the number of complete BUSCOs is 2,383, or 97.5%. We chose to remove contigs smaller than 750 bp in all of our assemblies (Table 5) because the combined read length for forward and reverse reads is 500 bp and anything smaller than 750 bp is not likely to be of high quality.

The assessment of the contamination of the P. coquilletti de novo genome using Blobology (Kumar et al., 2013) reveals that there is very little contamination with the vast majority of hits being arthropods (57.46%) or either unmapped (21.04%) and no-hit (21.19%) (Figs. S1–S2 and Table S2).

Preliminary annotation with MAKER (Drosophila reference libraries) produced 10,246 putative genes. D. melanogaster has a comparable genome size (164 Mbp compared to 210 Mbp for P. coquilletti; NCBI), but more than 17,000 genes (FlyBase.org). We plan to refine our preliminary annotation with more training and manual curation in order to improve our estimate. The P. coquilletti MAKER GFF file has been deposited at Figshare (doi: 10.6084/m9.figshare.4055643).

After we finished analyzing both the transcriptomes and genome, and conducted the phylogenetic analyses discussed below, we made a tally of the software programs used to generate these results: 30. This does not include all dependencies, so it is a bit of a conservative estimate. The sheer number of pieces of software in which a researcher using genome-scale data must feel at least conversant is quite large. These data do not lend themselves well to bioinformatics pipelines, either, because there are constantly improvements in existing software that change something about their usage or even new software that is found to be better for certain portions of the workflow and flexibility is key, which can come at the expense of a fairly steep learning curve for researchers just getting their feet wet with genome-scale data. The best way to combat the barrier (whether it is perceived or real) for those who might be interested in developing a set of genomic resources for non-model taxa is to thoroughly document which and how existing software was used (Fig. 2). While it is not possible to exactly recreate one group’s analysis because the computational infrastructure available is undoubtedly different, keeping track of how analyses are done is the first step toward reproducible work. Finally, we know that because tools are constantly improving, all current genomic resources are really just drafts that will be improved upon as we improve the software and databases upon which our results rely.

Phylogenetic trees

OMA produced 9,080 putative orthologs shared among four or more species. There were only a small number of loci found in all taxa because of the four Ion Torrent transcriptomes are of subpar quality. While we plan to resequence these on an Illumina platform in the future, they did produce enough data to be placed convincingly in a phylogenetic framework, which will be used to design exon capture probes in order to sample hundreds more species for which we have DNA but no RNA-quality specimens. The concatenated alignment length was 4,936,984 amino acid residues. Individual locus alignments as well as the concatenated alignment and best partitioning scheme from PartitionFinderProtein have been deposited at Figshare (concatenated alignment doi: 10.6084/m9.figshare.4055622, PartitionFinderProtein doi: 10.6084/m9.figshare.4055814).

Orthologous loci obtained from the transcriptomes were used to construct phylogenetic trees (Figs. 5–6) for the 16 taxa. While the small taxon sampling is not sufficient to provide new insights into the relationships within Mydidae (3 species included) or Asilidae (10 species), some comments on the higher-level relationships among and within families can be made. Both hypotheses (Figs. 5–6) support the position of Bombyliidae (Thevenetimyia californica) being more closely related to the other Asiloidea taxa Apioceridae, Asilidae, and Mydidae than to Tabanidae (Tabanomorpha, see also Fig. 1). A taxon Apioceridae plus Mydidae is supported as monophyletic and as the sister-group to Asilidae as previously proposed (Dikow, 2009b; Dikow, 2009a; Trautwein, Wiegmann & Yeates, 2010; Wiegmann et al., 2011). The relationships within Mydidae support the monophyly of the subfamily Mydinae with the two included genera Messiasia and Mydas. Within Asilidae, the four included Asilinae genera form a monophylum (Machimus, Philonicus, Proctacanthus, and Tolmerus) with Proctacanthus recovered as sister-group to the remaining three genera. The clade (Laphystia (Diogmites + Nicocles)) is recovered in both analyses, but with a different sister-group. The placement of the genus Lasiopogon, which has the fewest data available (Table 3), differs in the analyses. The tree based on concatenated loci (Fig. 5) places Lasiopogon as sister-group to Scleropogon, which is the least supported clade (bootstrap value of 36), and both genera are placed together as sister-group to (Laphystia (Diogmites + Nicocles)). In the ASTRAL analysis (Fig. 6), Lasiopogon is placed as sister-group to the Asilinae. The Laphriinae, here represented by Laphystia, has not been recovered as the earliest divergence within Asilidae has postulated by Dikow (2009b) and Dikow (2009a).

Time-calibrated tree

†Araripogon axelrodi and †Cretomydas santanensis, representing the oldest, definitive Asilidae (Grimaldi, 1990) and Mydidae (Willkommen & Grimaldi, 2007) fossils, respectively, were placed as sister group to the remaining taxa from these families. †Burmapogon bruckschi, a 100 myo Burmese amber fossil (Dikow & Grimaldi, 2014), was placed as sister-group to Lasiopogon (see discussion of placement in Dikow & Grimaldi (2014)). Four species of Asilidae are known from Baltic amber with an age of 45 my (Evenhuis, 1994) and a study of 25 newly discovered amber specimens (T Dikow, 2015, unpublished data) adds three additional species. In total, four Baltic amber assassin flies (two previously described and two yet undescribed) are included as they can be sufficiently well-placed. Based on the most recent phylogeny of Asilidae (Dikow, 2009b), †Asilus klebsi was placed as sister-group to Tolmerus + Machimus + Philonicus, a yet undescribed Asilinae genus and species with an ovipositor similar to that of extant Promachus (Asilinae: Apocleini) as sister-group to the included four Asilinae species, the Laphriinae: Atomosiini †Protoloewinella keilbachi as sister-group to Laphystia tolandi, and an undescribed genus and species of Stenopogoninae as sister group to Scleropogon duncani.

With the inclusion of seven Cretaceous and Tertiary fossils, we provide the first time-calibrated tree for ages within Apioceridae, Asilidae, and Mydidae (Fig. 7). The limited taxon sampling in our analysis prohibits a detailed discussion of clade ages, nonetheless the results provide a first view of the earliest divergences within these taxa. The split between Apioceridae + Mydidae and Asilidae is postulated to have occurred about 158 Million years ago (mya), which is approximately 25 my earlier than hypothesized by Wiegmann et al. (2011). Likewise, the split between Apioceridae and Mydidae is here postulated to be 135 mya, 15 my earlier than the age estimated by Wiegmann et al. (2011). Within the included Mydidae, the split between Ectyphinae and Mydinae is postulated to have occurred 107 mya, which is supported by the placement of †Cretomydas (112 myo) before the Ectyphinae/Mydinae divergence in a morphological phylogenetic study on the family (T Dikow, 2015, unpublished data). Based on our analysis, the earliest divergence within Asilidae occurred 128 mya and that of the four included Asilinae genera 81 mya. The remaining five subfamilies diverged 114 mya, which is supported by the placement of †Burmapogon (100 myo) within this clade.