A scaffold-level genome assembly of a minute pirate bug, Orius laevigatus (Hemiptera: Anthocoridae), and a comparative analysis of insecticide resistance-related gene families with hemipteran crop pests

Background Orius laevigatus, a minute pirate bug, is a highly effective beneficial predator of crop pests including aphids, spider mites and thrips in integrated pest management (IPM) programmes. No genomic information is currently available for O. laevigatus, as is the case for the majority of beneficial predators which feed on crop pests. In contrast, genomic information for crop pests is far more readily available. The lack of publicly available genomes for beneficial predators to date has limited our ability to perform comparative analyses of genes encoding potential insecticide resistance mechanisms between crop pests and their predators. These mechanisms include several gene/protein families including cytochrome P450s (P450s), ATP binding cassette transporters (ABCs), glutathione S-transferases (GSTs), UDP-glucosyltransferases (UGTs) and carboxyl/cholinesterases (CCEs). Methods and findings In this study, a high-quality scaffold level de novo genome assembly for O. laevigatus has been generated using a hybrid approach with PacBio long-read and Illumina short-read data. The final assembly achieved a scaffold N50 of 125,649 bp and a total genome size of 150.98 Mb. The genome assembly achieved a level of completeness of 93.6% using a set of 1658 core insect genes present as full-length genes. Genome annotation identified 15,102 protein-coding genes - 87% of which were assigned a putative function. Comparative analyses revealed gene expansions of sigma class GSTs and CYP3 P450s. Conversely the UGT gene family showed limited expansion. Differences were seen in the distributions of resistance-associated gene families at the subfamily level between O. laevigatus and some of its targeted crop pests. A target site mutation in ryanodine receptors (I4790M, PxRyR) which has strong links to diamide resistance in crop pests and had previously only been identified in lepidopteran species was found to also be present in hemipteran species, including O. laevigatus. Conclusion and significance This assembly is the first published genome for the Anthocoridae family and will serve as a useful resource for further research into target-site selectivity issues and potential resistance mechanisms in beneficial predators. Furthermore, the expansion of gene families often linked to insecticide resistance may be an indicator of the capacity of this predator to detoxify selective insecticides. These findings could be exploited by targeted pesticide screens and functional studies to increase effectiveness of IPM strategies, which aim to increase crop yields by sustainably, environmentally-friendly and effectively control pests without impacting beneficial predator populations. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08249-y.


Background
Loss of crops to insect pests can account for ~ 10% of potential yield, as a result of both direct feeding damage and the transfer of viral plant diseases [1]. Thus, to maximise crop yields and sustain food production for a growing world population, pests need to be controlled. At present this control relies mainly on the use of synthetic pesticides, many of which are non-selective and are therefore toxic to both their target pest species and to beneficial predators and parasitoids. As a result there may be a reduction in the predator populations to a level where they are no longer able to contribute natural pest control. This, along with the development of insecticide resistance in pests, can lead to pest populations surging, sometimes to even higher levels than pre-pesticide application [2][3][4]. Beneficial predators, such as those in the genus Orius, have proven to be especially effective in the biological control of crop pests [5]. As generalist predators, Orius species target a wide variety of pest species including aphids, beet armyworm, leafhoppers, mites, thrips and whiteflies, many of which are the world's most damaging crop pests [6,7]. Some Orius species are commercially available as biological control agents and are widely used for this purpose as part of integrated pest management (IPM) strategies, especially in covered crops [8][9][10].
Whole genome sequences of insects are helping us to understand many aspects of their biology and behaviour, and this can be applied to potential insecticide resistance mechanisms in pest insects and their natural enemies. However, only a few genomes of beneficial predator species have been published to date, including a phytoseiid mite, Galendromus occidentalis [11]; three parasitoid wasps, Nasonia giraulti, Nasonia longicornis and Nasonia vitripennis [12] and two lady beetles, Harmonia axyridis and Coccinella septempunctata [13]. To date there are no published genomes for species of the Hemiptera: Anthocoridae (i.e. minute pirate bug) family of predators. In contrast, a growing number of genomes of crop pests are available [14][15][16][17][18][19][20][21][22][23][24][25][26]. This larger number of pest genomes, relative to beneficial predator genomes could be in part because up until recently, the genomes of the pests themselves have appeared more useful in terms of developing targeted pesticides and investigating mechanisms of pesticide resistance. However, agriculture is now moving increasingly away from pesticide use -particularly with the Directive on Sustainable Use of Pesticides 2009/128/EC [27] -and towards IPM strategies, which includes the use of beneficial predators. Future studies of pesticide resistance mechanisms should therefore include beneficial predator genomes alongside pest genomes in order to help select targeted pesticides which do not harm beneficials and subsequently improve the efficacy of IPM strategies [28][29][30][31][32].
The aim of the work reported here was to develop a high-quality genome assembly for O. laevigatus, to serve as a resource for research into this species as well as the wider Anthocoridae family, which consists of 400-600 mostly predaceous minute pirate bug species -a potentially valuable source of biological control agents [33]. The O. laevigatus genome was then used for comparative analyses between beneficial predators and crop pests, focusing on genes encoding potential insecticide resistance mechanisms.
There are two main types of insecticide resistance mechanisms: increased expression of genes encoding protein families involved in metabolic resistance and point mutations in genes encoding insecticide target proteins [34]. Gene families involved in insecticide resistance in pest species are known to include cytochrome P450 monooxygenases (P450s), ATP binding cassette transporters (ABCs), glutathione S-transferases (GSTs), UDP-glucosyltransferases (UGTs) and carboxyl/cholinesterases (CCEs) [35][36][37][38][39][40]. Comparisons of the genes/ proteins which may be involved in insecticide resistance in crop pests with the corresponding genes in beneficial insects, could aid the development of insecticides which target crop pests but have limited impact on beneficial predator populations. This could prove key to developing successful IPM strategies which exploit differences in the ability of predators and crop pests to tolerate pesticides. Improving the availability of beneficial predator genomes could also help the selection of beneficial predators with genes/mutations for inherent insecticide resistance before being released in the field for biological control [41].
The results presented here provide a comprehensive foundation for further study of potential insecticide resistance mechanisms in beneficial predators and how they compare to crop pests.

Sample preparation and sequencing
Orius laevigatus (commonly known as a minute pirate bug) were obtained from 'Bioline AgroSciences' . CO 2 was used for anaesthesia to allow the insects to be sorted from the substrate. Both adults and nymphs were then flash frozen with liquid N 2 and stored at − 80 °C. The whole process was done within 48 h of arrival. 1000 individuals were pooled for genomic DNA/ RNA extractions, which were carried out in-house at Rothamsted Research. The commercial DNAzol reagent was used for the DNA extractions, and the Bioline Isolate II RNA Mini Kit was used for the RNA extractions. The DNA and RNA were sent for library preparation and sequencing by Genewiz (New Jersey, US).
The genome assembly was developed using a hybrid assembly strategy with both Illumina short reads and Pacific Biosciences (PacBio) long reads.
Short reads were sequenced using 2 mg of DNA and a library with an insert size of 200 bp. Sequencing was done using Illumina HiSeq 4000 with a 2x150bp pairedend configuration. 413,143,574 reads were obtained with a total length of 123 Gb (820x). Raw reads are available under SRA accession: ERR6994870. K-mer counting of the raw Illumina DNA data was done using Jellyfish 2.2.6 [42]. Canonical (−C) 21-mers (−m 21) were counted and a histogram of k-mer frequencies outputted. Genom-eScope 2.0 [43] was used to process this histogram with 'ploidy' set at 2 and 'maximum k-mer coverage cut-off ' set at 10,000.
To obtain long read PacBio data, 3.7 mg of DNA first underwent blue pippin size selection (> = 10 kb) to remove low molecular weight DNA. < 500 ng of DNA remained after size selection, and so a low input protocol was used for library construction with an insert size of 20 kb. Sequencing was done using the PacBio Sequel I platform and 537,651 reads were obtained with a total length of 6Gb (44x) and an N50 of 11,287 bp. Raw reads are available under SRA accession: ERR6941611.
Transcriptome sequencing used 10 mg of RNA and a library construction with an insert size of 150 bp and PolyA selection for rRNA removal. Sequencing was done using Illumina HiSeq 4000 with a 2x150bp paired-end configuration. 413,137,378 reads were obtained. Raw reads are available under SRA accession: ERR7012629.
FastQC v.0.11.8 [44]. was used for quality checks on the raw Illumina HiSeq DNA and RNA sequence data. Adapters were trimmed, low-quality bases (below a score of 3) were removed from the start and end of reads and any reads with a length less than 36 bases were also removed. Trimmomatic v.0.38 [45]. was used for these trimming steps. Quality trimming of reads using Trimmomatic resulted in a 0.2% loss of reads for whole genome sequencing and a 5% loss of reads for transcriptome sequencing (Table 1).

Genome quality assessment
Basic metrics from the genome assembly were calculated using a script developed for the ' Assemblathon' [46]. These metrics include scaffold/contig N50, longest and shortest scaffold length, number of scaffolds exceeding a range of lengths and number of gaps/N's in the assembly.
The completeness of the genome assembly and annotation for Orius laevigatus was assessed using the Benchmarking Universal Single-Copy Orthologs (BUSCO) [47] of the insect gene set (insecta odb9). 'Genome' mode was used to assess the assembly, and 'protein' mode to assess the annotation. 'Fly' was used as the training species for Augustus gene prediction. BUSCO assessments were then run with default parameters.

De novo genome assembly
The overall assembly pipeline is shown in Fig. 1. The raw PacBio long reads were assembled into contigs with the Flye v2.5. de novo assembler [48,49]. Rascaf was then used to improve the Flye genome assembly with RNAseq data [50]. Contigs were also produced with the raw PacBio long reads using Canu v1.8 [51] as well as with FALCON v1.3.0 and FALCON-Unzip, which is recommended for heterozygous/outbred organisms with diploid or higher ploidy (and also includes phased-polishing with Arrow) [52,53].
QuickMerge v0.3 [54] was used to merge the assemblies, with Flye as the reference assembly. BUSCO outputs were compared between the merged assembly and the standalone assemblies to identify genes which had been lost during the merging process. Full-length contigs containing these missing genes were extracted from the standalone assemblies and added to the merged assembly, based on the assumption that these contigs would also contain other missed genes (i.e. those not included in BUSCO's list of 1658 core insect genes). Multiple rounds of Pilon error polishing [55] were performed, using the Illumina short read data, until no further improvement in BUSCO score was seen. Redundans [56] was used for scaffolding and redundant contig removal. Redundans is geared towards highly heterozygous genomes. Some redundant regions had to be removed manually, as Redundans does not detect redundancy when only part of the contig is duplicated. The nucmer tool from the MUMmer4 package [57] was used to detect these redundant regions through a whole genome self-alignment.
A BLAST search against the NCBI Reference Sequence (Refseq) database release 93 [58], was performed using the Tera-BLAST algorithm on a Time-Logic DeCypher system (Active Motif Inc., Carlsbad, CA). The results were processed with Megan [59] to identify any bacterial or viral sequences which were then removed manually in Geneious v10.2.6.
The mitochondrial genome sequence was identified and extracted by running a BLAST search of the O. laevigatus genome against the Orius sauteri mitochondrial genome which is publicly available at NCBI, Gen-Bank accession No. KJ671626 [60].
A de novo species specific repeat library was constructed using RepeatModeler v1.0.7 [66] to identify repeat models. These models were searched against the GenBank non-redundant (nr) protein database for Arthropoda (e value < 10 − 3 ) using Blastx to remove any potential protein-coding genes. This was combined with transposon data to create a custom library. Transposons were identified from the transcriptome assembly by running HMMER: hmmscan [67] against the Pfam database [68] and filtering the resultant Pfam descriptions for those containing "transposon". A search for transposons was also done on transcripts produced from MAKER and these transposons were then added to the custom repeat Fig. 1 The assembly pipeline for the Orius laevigatus genome library which was used for a second round of MAKER. RepeatMasker v4.0.7 [69] was used to mask repeats in the genome assembly using these repeat libraries, as well as to estimate the abundances of all predicted repeats.
Evidence from assembled transcripts was transferred to the genome assembly via MAKER. The output from this was then used to produce a high confidence level gene model training set -overlapping and redundant gene models were removed. Augustus and GeneMark were trained using this training set prior to being used for ab initio gene predictions. FGeneSH was run based on the Drosophila melanogaster genome.
The best transcripts from both the ab initio gene prediction annotation and the transcriptome-based annotation were selected using Evigene (classified by reasonable protein size and homology to other species) and combined to create the final annotation.

Comparative genomics and phylogenetic analysis
To produce the species tree, orthogroup gene trees were produced using Orthofinder [77] and the tree was inferred from these using the STAG method [78].
In order to identify genes potentially involved in insecticide resistance, the PFAM domains assigned to gene models during annotation (as described in the 'Genome Annotation' methods section) were used as follows: CCEs (PF00135/IPR002018), GSTs (IPR004045/PF02798), (IPR004046/PF00043), P450s (IPR001128/PF00067), ABCs (IPR003439/PF00005) and UGTs (IPR002213/ PF00201). Proteins from UniProt for the classes of interest, from hemipteran species, were used for BLAST queries against O. laevigatus to identify any missed genes and to assist with subfamily assignment within these classes. Subfamily assignment for O. laevigatus gene families was finalised using phylogenetic trees produced using MAFFT alignments [79,80] and RaxML v8.2.11 [81]. The GAMMA LG protein model [82] was used and a bootstrap consensus tree was inferred from 100 replicates.
Manual checks and curation were performed for genes potentially involved in insecticide resistance. Increased copy numbers of these genes often led to adjacent tandem duplications being incorrectly annotated as one gene model, therefore curation was important to prevent incorrect gene numbers being reported in later analyses. The exon/intron boundaries and start/stop codons of the genes were confirmed through visualization in IGV [83] of RNAseq data mapped to the genome using HISAT2 v2.0.5 [70] and the gene models were edited in Geneious where necessary.
The P450s were classified and named by Dr. David Nelson [84]. The UGTs were classified and named by Dr. Michael Court [85]. Nomenclature of P450s and UGTs is based on the evolutionary relationships of the sequences. P450 and UGT sequences were BLAST searched against named insect sequences and were assigned to known families if they were > 40% (for P450 families) or > 45% (for UGT families) identical. Other sequences were assigned to new families based on their clustering on trees and their percent identity to each other.

Sequencing
In order to produce enough DNA and RNA for sequencing, ~ 1000 individuals of O. laevigatus were required. Because they were obtained commercially, the level of inbreeding of the culture was not known. However, all individuals were obtained from a single colony within the rearing facility. A high heterozygosity level was therefore a possibility and this was kept in mind when making decisions during the assembly process.

Genome metrics evaluation based on raw reads
The raw read k-mer analysis with GenomeScope 2.0 estimated a haploid genome size of ~ 141 Mb (Table 2), in line with the final assembly size of 151 Mb. A genome size estimate using methods such as flow cytometry would have provided a more accurate estimate, however, such data was not available for the Orius genus. This could be considered a limitation to the study, as 141 Mb was provided as a genome size estimate to Canu, Flye and FALCON-Unzip which may have affected the outputted assemblies. Genome repeat length was 20 Mb, 16.5% of the total estimated genome size.
The heterozygosity rate ranged from 1.20 to 1.30%. This alongside the small 'shoulder' to the left of the main 'full-model' peak (Fig. 2), indicates a fairly high level of heterozygosity, which was taken into consideration in the assembly strategy.

Assembly
Flye, FALCON and Canu were used to produce 3 separate genome assemblies. The statistics for these assemblies, as well as for subsequent versions of the assembly outlined in this section are shown in Additional file 2. Rascaf improved the contiguity of the Flye assembly through alignment of the RNA-seq data to the genome, likely because it is less affected by the use of multiple individuals versus genome assembly tools which include non-conserved sequences from a population of individuals. FALCON-Unzip improved the FALCON assembly contiguity with a 4.5-fold decrease in the total number of scaffolds (although this coincided with a ~ 9% loss of complete gene models found using BUSCO and suggests that FALCON-Unzip may have been too stringent for this genome -perhaps because it was designed with plant and fungal genomes in mind [52,53]).
Flye (both with and without Rascaf ) had the best assembly statistics in terms of scaffold N50 and BUSCO score. However, FALCON-Unzip achieved the largest 'longest scaffold' of the three assemblers.
Quickmerge was used to merge the FALCON-Unzip assembly, Rascaf improved Flye assembly and the Canu assembly. The resultant merged assembly had better continuity than any of the stand-alone assemblies, however, the BUSCO completeness was slightly worse than the standalone Flye assembly (and worsened with the second round of Quickmerge). This was likely due to mis-assemblies in the component assemblies causing alignment issues, which resulted in sections of the misassembled contigs being discarded.
Pilon was used for error polishing and improved the BUSCO completeness score. Redundans (redundancy removal and scaffolding/gap-closing) improved the scaffold N50 and removed redundant scaffolds.
A comparison of the gene models (core insect genes from the insecta odb9 BUSCO gene set) found in the original Flye / FALCON-Unzip / Canu assemblies versus the merged assembly showed that some of the gene models were found in at least one of the original assemblies, but were missing in the merged assembly. Of the 154 missing or fragmented genes in the merged assembly (out of a total 1658 core insect genes), 5 were found in the FALCON-Unzip assembly, 5 in the Flye assembly and 46 in the Canu assembly. Manual editing to bring the full-length contigs containing these missing genes into the merged assembly took the BUSCO completeness score up by 5%. A final round of Pilon improved this score by an additional 0.5% (further rounds of Pilon did not improve the score).

Annotation
Gene prediction with MAKER identified 15,102 proteincoding genes with the encoded proteins having a mean length of 464 amino acids. Of these, 12,949 (86%) had a match to NCBI's non-redundant (nr) database and 11,616 (77%) contained InterPro motifs, domains or signatures. In total, 13,112 (87%) were annotated with either blastp or InterPro and 10,192 were annotated with a GO ID. More information on the InterPro member database annotations is given in additional file 1. The longest protein found was an 'egf-like protein' at 14,628 amino acids. The resultant gene set was 84.5% BUSCO (insecta) complete.
From the Infernal tool inference of RNA alignments, a total of 791 non-coding RNA elements and 269 cis-regulatory elements were found in the genome (Table 4).

Repeat annotation
Transposable and repetitive elements made up 27.07% of the assembled O. laevigatus genome (Table 5) and the majority of these (20.4%) were unclassified repeats. This is close to the reported repeat content of other hemipteran species, for example: Cimex lectularius -31.63% [86] and Acyrthosiphon pisum -38% [15], an exception is Rhodnius prolixus which has an unusually low repeat content of 5.6% [87].

Mitochondrial genome
A circularized mitochondrial genome of 16,246 bp, assembled and annotated using MITOS2, consisted of 13 protein coding genes, 19 tRNA genes, 2 rRNA genes and an A + T rich region with a length of 1460 bp and  an A + T content of 72.7% (Fig. 3). This closely matches the Orius sauteri mitochondrial genome, which is also 16,246 bp and has 13 protein-coding genes, 22 tRNA genes, 2rRNA genes and an A + T rich region of 1758 bp and an A + T content of 73.5% [60].

Phylogeny
OrthoFinder assigned 318,985 genes (88.8% of total) to 27,481 orthogroups. There were 1621 orthogroups with all species present and 45 of these consisted entirely of single-copy genes.
Phylogenetic analysis correctly clustered O. laevigatus within the hemipteran clade (Fig. 4) and identified Cimex lectularius as its closest relative.

Comparative genomics ABC transporters
ATP-binding cassette transporters (ABCs), the largest known group of active transporters, can eliminate xenobiotic compounds -such as secondary metabolites produced by plants or insecticides -through translocation [36]. These transporters are subdivided into eight subfamilies: ABCA-H. ABCB, ABCC and ABCG are the subfamilies most associated with resistance to a variety of insecticides including pyrethroids, carbamates, organophosphates and neonicotinoids [88]. Forty-one of the 64 transporters in O. laevigatus belong to these 3 class-specific expansions ( Table 6) which could confer resistance to insecticides (a phylogenetic tree showing relationships of ABC transporters in O. laevigatus is included in Additional file 3, full sequences are included in Additional file 4). Table 6 shows a comparison of numbers of ABC transporter genes found in the current study with those reported for some pest species. The gene family expansions were generally seen in the ABCC and ABCG classes for all hemipteran species and slightly larger expansions were seen in some crop pests compared to O. laevigatus for the ABCC class, however, the expansions were of very similar size for both crop pests and O. laevigatus in the ABCG class. Overall, the total numbers of ABC transporter genes were similar across all the hemipteran species compared.

Glutathione S-Transferases
The glutathione S-transferases (GSTs) protein family is large and functionally diverse, and is known to confer resistance to all main insecticide classes. GST-mediated detoxification of insecticides takes place via several different mechanisms, including protecting against   oxidative stress, binding and sequestration of the insecticide, and by catalysing the conjugation of glutathione to the insecticide to reduce their toxicity [37]. The number of GST genes in O. laevigatus was fairly similar to other hemipteran close relatives, with the exception of the sigma class, which was notably larger ( Table 7, full sequences included in Additional file 4). Of the 16 genes in the sigma class, 9 genes (mRNA13082 and mRNA13086-13,093) were adjacent on the same scaffold, indicating a lineage specific expansion (Fig. 5). Expansions in this class have been reported in several hemipteran species including Triatoma infestans, Myzus persicae, Halyomorpha halys and Murgantia histrionica [23,24,95,96]. The sigma class has been found to play an important role in detoxification of organophosphorus insecticides in hemipteran species [97], therefore this expansion could potentially confer some tolerance to organophosphates in O. laevigatus. The delta and epsilon classes of GSTs are linked to insecticide resistance to pyrethroids [98,99]. The delta class is much larger in several crop pests compared to O. laevigatus and its close relatives which suggests these crop pests could exhibit a higher level of delta class GST-mediated pyrethroid resistance. The epsilon class has previously been thought to be specific to Holometabola [100], and whilst Trialeurodes vaporariorum has a single member, the epsilon class is absent from all other Hemiptera species, suggesting potential epsilon class GST-mediated pyrethroid resistance is most likely absent in O. laevigatus and its close relatives, as well as most Hemiptera crop pests.

Carboxyl/cholinesterases
Many carboxyl/cholinesterases (CCEs) are linked to detoxification of organophosphorus, carbamate and pyrethroid insecticides and acetylcholinesterase (AChE) is the target for organophosphate and carbamate insecticides, with amino acid substitutions being linked to resistance [39]. Thirty-two members of the CCE superfamily, including 1 AChE gene, were found in the O. laevigatus genome (Table 8, full sequences included in Additional file 4) which is a similar number to that reported for Cimex lectularius, which had 30 CCE genes and 2 AChE genes [86].
The dietary class of CCEs is involved in insecticide and xenobiotic detoxification [106]. O. laevigatus has no genes within this class, in line with T. infestans and   The dietary class is involved in pyrethroid resistance [107]; however, T. infestans exhibits pyrethroid esterase activity despite having no dietary esterases [108]. O. laevigatus has also shown the ability to develop pyrethroid resistance -although the exact mechanism of this resistance is not yet known [109]. The hormone and semiochemical processing class is also involved in insecticide metabolism, due to the presence of β-esterases [110,111]. There may be some redundancy in genes potentially involved in insecticide detoxification from the dietary and hormone/semiochemical processing classes. This might explain why only one of these classes shows an increased number of genes for each of these hemipteran species (Table 8), as having increased numbers of both classes would be redundant, whilst very low numbers of both classes would be detrimental. The lack of the dietary class may therefore not impact the xenobiotic resistance abilities of O. laevigatus, as it has 16 genes within the hormone/semiochemical processing class.
The remaining CCEs in O. laevigatus belong to the neurodevelopmental class and include the neuroligins, gliotactins, glutactins and neurotactins, which are noncatalytic due to the lack of a critical serine residue. Acetylcholinesterase is the only protein in this class which has been linked to organophosphate resistance [112,113].

UGTs
UDP-glucosyltransferases (UGTs) are detoxification enzymes speculated to be involved in insecticide metabolism. Although the exact mechanisms of UGT-mediated resistance have not yet been identified, their upregulation has been shown in resistant strains of P. xylostella [35] and they have been linked to diamide resistance in Chilo suppressalis [114], neonicotinoid resistance in Diaphorina citri [115] and they also contribute to insecticide detoxification via the elimination of oxidative stress in Apis cerana [116].
The number of UGT genes in O. laevigatus was much lower than for other hemipteran species (Table 9). The UGTs were submitted to Dr. Michael Court for naming and are included in Additional file 4. Numbers of UGTs have been reported to be lower in non-phytophagous  Table 8 Numbers of Carboxyl/cholinesterases annotated in Orius laevigatus (this study), Cimex lectularius [101], Rhodnius prolixus [103], Triatoma infestans [96], Frankliniella occidentalis [90], Myzus persicae [95], Acyrthosiphon pisum, Bemisia tabaci [104] and Trialeurodes vaporariorum [105] and their distribution across classes and clades * C. lectularius numbers may be an underestimate as sequencing coverage was low for this study, clade assignment was also uncertain as a result ** A more recent study [86] found 30 CCE genes in C. lectularius, and is more likely to be a true representation, but they had not been assigned into classes/clades *** Numbers in brackets represent the possible true numbers of R. prolixus CCEs, based on a potential misassignment insects [91], which could explain the low numbers seen in O. laevigatus and R. prolixus compared to crop pests. This suggests that UGT-mediated detoxification may be lower in O. laevigatus than in crop pests.

Cytochrome P450s
Cytochrome P450s are a diverse superfamily capable of metabolizing a huge variety of endogenous and exogenous substrates. In insects they are associated with growth and development, metabolism of pesticides and plant toxins as well as the production and metabolism of insect hormones and pheromones. P450s are associated with resistance to insecticides from a variety of classes, including pyrethroids, carbamates and neonicotinoids. They are also linked to the activation of organophosphates and other pro-insecticides (a form of insecticide which is metabolized into an active form inside the host) [38]. Upregulation of P450s in insects has been shown to confer insecticide resistance [119][120][121][122], and conversely downregulation occurs in response to pro-insecticides [123,124]. A total of 58 full-length P450 genes were identified in the O. laevigatus genome, 11 P450 fragment genes were also found as well as 1 pseudogene. (Sequences are included in Additional file 4). These sequences were named by Dr. David Nelson using his in-house pipeline [84]. The majority of these genes (34) belonged to the diverse CYP3 class, which was a similar size to other hemipteran species (Table 10).
The CYP3 clade is currently the P450 clade most associated with insecticide resistance -notably the CYP6 and CYP9 families [128]. Interestingly the CYP9 family was not present in O. laevigatus, as found for T. infestans, R. prolixus, M. histrionica and H. halys [23,24,96]. Further investigation into the assignment of classes within the CYP3 clade suggests the lack of the CYP9 class could be a common feature within Hemiptera (Table 10).
Expansion of the CYP397 gene family was seen in O. laevigatus, (Fig. 6) with 7 full-length CYP397 genes and 1 fragment. CYP397B1, CYP397B2, CYP397B6 and CYP397C1 were directly adjacent on the same scaffold, indicating tandem duplications. Sequence similarity of the CYP397 genes to CYP397B1 ranged from 52 to 86%, which suggests a variation in ages of these tandem duplications. Cimex lectularius also showed an increased copy number of CYP397 with 6 copies (A1-A6) [86]. CYP397A1 is significantly upregulated (> 36 fold) in pyrethroid-resistant strains of C. lectularius [127], therefore the expansion of this gene family could potentially confer some tolerance to pyrethroids in O. laevigatus.
A previous study [129] looked at the effect of insecticide synergists on Orius tristicolor (another minute pirate bug of the Anthocoridae family), and found that PBO (an inhibitor of P450s and esterases) significantly increased the mortality rate when combined with indoxacarb (an oxadiazine insecticide). Whereas inhibition of solely GSTs or esterases did not reduce mortality. Upregulation of P450s, esterases and GSTs have all been seen in response to oxadiazines [130], therefore the fact that only P450 inhibition had an impact on mortality rate suggests P450s may be the primary detoxification mechanism of O. laevigatus.

Target site mutations
Point mutations resulting in amino acid substitutions in the target proteins of insecticides have been characterised in many insecticide resistant insect species, including in the sodium channel gene para which confers resistance to pyrethroids [131]; the acetylcholinesterase-1 (ace-1) enzyme associated with organophosphate resistance [132] and the acetyl-coenzyme A Carboxylase (ACC) enzyme linked to keto-enol (spirotetramat) resistance [133]. Despite these mutations having been observed in a variety of hemipteran crop pests, none were observed in this O. laevigatus assembly. Although, it is important to note that the O. laevigatus assembly was a consensus of ~ 1000 individuals, therefore differences in target sites would likely only be apparent if they were present in the majority of the population. Overall, tolerance of insecticides by O. laevigatus resulting from target site differences seems unlikely compared to what is seen in crop pests, where there has been intensive selection pressure.
The ryanodine receptor (RyR) is the target of diamide insecticides, and two target site resistance mutations conferring amino acid substitutions (I4790M and G4946E -numbering according to Plutella xylostella, PxRyR) have been identified in lepidopteran pests [21,134]. Interestingly, O. laevigatus has the I4790M substitution (full sequence for Ryanodine in Additional file 4) which has been shown to confer varying levels of resistance to diamides. This point mutation was Table 9 Numbers of UDP glucuronosyltransferase genes found in O. laevigatus (this study), Rhodnius prolixus, Tetranychus urticae, Nilaparvata lugens, Acyrthosiphon pisum, Bemisia tabaci [19], Myzus persicae [117] and Trialeurodes vaporariorum [118] Table 10 Total numbers of Cytochrome P450 genes annotated in Orius laevigatus (this study), Cimex lectularius [86], Rhodnius prolixus, Triatoma infestans [96], Frankliniella occidentalis, Thrips palmi [90], Myzus persicae, Acyrthosiphon pisum [95], Trialeurodes vaporariorum [93], Bemisia tabaci [125], Halyomorpha halys [23] and Murgantia histrionica [24] *Values in brackets represent total gene numbers including partial and fragment genes. For other species partial and fragment p450 genes were excluded in cases where they were listed as such -some may remain in the counts if official naming and curation had not taken place **Values used are those from [86], but values differed by study [126]   also present in other hemipteran species as shown in Fig. 7 (except for Lygus hesperus which had an I > L mutation). I4790M has been detected in lepidopteran populations across the globe and is considered to be a 'selectivity switch' for diamides [135]. O. tristicolor showed high levels of resistance to chlorantraniliprole (a diamide insecticide) with < 5% mortality [129] with the I4790M substitution being the main cause [136]. It is therefore possible that I4790M may confer some tolerance to diamides in O. laevigatus, and indeed, diamide resistance has been reported in O. laevigatus [137]. However, I4790M could potentially also confer diamide tolerance in crop pests -diamide resistance has already been shown in F. occidentalis [137]. Therefore this would likely not be an exploitable difference for IPM strategies.

Conclusions
PacBio long-read technology combined with low errorrate short-read Illumina sequencing enabled the production of a high-quality genome and mitochondrial assembly for O. laevigatus. Whilst genome continuity may not be as good as an assembly generated from a single insect, the genome completeness is still of a sufficient quality to aid with comparative and functional genomics analyses and provides a useful first reference genome for the Anthocoridae family. An experimental estimate to confirm genome size and Hi-C based scaffolding would likely be the next best steps to significantly improve this genome in the future. Comparative analyses of O. laevigatus with hemipteran crop pests showed evidence of possible differences in xenobiotic tolerance, including a potential increase in GST-mediated tolerance of organophosphates in O. laevigatus, whilst GST-mediated pyrethroid tolerance may be more prevalent in crop pests. There may also be less UGT-mediated tolerance to diamides and neonicotinoids in O. laevigatus compared to crop pests -although, the I4790M target site mutation may confer some degree of diamide insensitivity to O. laevigatus.
A recent study shows that there is significant variation in the susceptibility of O. laevigatus to pyrethroids when Fig. 6 Phylogenetic tree of the Orius laevigatus cytochrome P450s. The Cyp397 gene family is a member of clan 3. Amino acid sequences were aligned using MAFFT and analysed using RAxML (the GAMMA LG protein model was used). The bootstrap consensus tree was inferred from 100 replicates a variety of wild and commercial populations are assessed [109]. This suggests that beneficial predators such as O. laevigatus are certainly capable of developing insecticide resistance, but a combination of factors result in resistance developing slower than in pest species. This could be due to beneficial predators having smaller population sizes, longer life cycles, less exposure to pesticides and a lack of continuous selection pressure -beneficial predators often need to be re-released each season as populations migrate to new areas in search of food sources. These differences will have resulted in a lesser degree of selection for resistance mechanisms in O. laevigatus and therefore any observed differences in potential sensitivity would only be at low levels. Further comparisons looking at differences in gene expansions, expression levels and key target site mutations between resistant and susceptible strains of O. laevigatus would provide more concrete evidence for the findings in this study.
In conclusion, this study indicates differences in potential mechanisms of resistance between crop pests and O. laevigatus which could be exploited when selecting targeted insecticides. An increase in the number of pesticides which are safe for beneficial predators such as O. laevigatus would be of significant impact to pest management, especially at a time when the list of pesticides effective against crop pests is growing ever shorter. The findings also suggest that O. laevigatus has the ability to develop resistance to a variety of insecticides which could be used to our advantage through the selective breeding and selection of resistant strains of O. laevigatus for use in IPM strategies.