Whole-genome sequencing analysis and protocol for RNA interference of the endoparasitoid wasp Asobara japonica.

Asobara japonica is an endoparasitic wasp that parasitizes Drosophila flies. It synthesizes various toxic components in the venom gland and injects them into host larvae during oviposition. To identify and characterize these toxic components for enabling parasitism, we performed the whole-genome sequencing (WGS) and devised a protocol for RNA interference (RNAi) with A. japonica. Because it has a parthenogenetic lineage due to Wolbachia infection, we generated a clonal strain from a single wasp to obtain highly homogenous genomic DNA. The WGS analysis revealed that the estimated genome size was 322 Mb with a heterozygosity of 0.132%. We also performed RNA-seq analyses for gene annotation. Based on the qualified WGS platform, we cloned ebony-Aj, which encodes the enzyme N-β-alanyl dopamine synthetase, which is involved in melanin production. The microinjection of double-stranded RNA (dsRNA) targeting ebony-Aj led to body colour changes in adult wasps, phenocopying ebony-Dm mutants. Furthermore, we identified putative venom genes as a target of RNAi, confirming that dsRNA injection-based RNAi specifically suppressed the expression of the target gene in wasp adults. Taken together, our results provide a powerful genetic toolkit for studying the molecular mechanisms of parasitism.

1. Introduction

Parasitism is an ecological lifestyle whereby the parasite organism exploits resources from the host organism to sustain life. In insects, parasitoid wasps belong to the order Hymenoptera, exploiting other insect species and arthropods as hosts. It is estimated that a group of parasitoid wasps accounts for more than 20% of all insect species. A large number of species have highly diversified strategies for successful parasitism on earth. For example, parasitoid wasps lay eggs at a certain developmental stage of their host arthropods, including the embryo, larva, pupa, and adult stages. Even in close relatives, some wasp species are adapted to many host species (generalist), whereas others only parasitize a limited number of host species (specialist). Many studies have reported ecological relationships and coevolution mechanisms between parasitoid wasps and their hosts. Previous studies have also revealed a wide variety of bioactive molecules, venom proteins, and virus-like particles produced in parasitoid wasps. Conversely, host arthropods acquire several repellant behaviours and immune defense systems against parasites., However, the underlying molecular mechanism of wasp parasitism is still poorly understood compared with the immune defense system of their hosts.

To study the molecular mechanisms of parasitism, it is crucial to obtain whole-genome information for both parasitoid wasps and their host arthropods. In the fruit fly Drosophila melanogaster, an excellent model organism, ∼50 drosophilid parasitoids have been reported., Among them, the genus Leptopilina (Hymenoptera: Figitidae) has been the best studied. These wasps lay their eggs in fruit fly larvae and the developing wasp larvae eventually consume and kill their fly hosts. Leptopilinaheterotoma is a generalist, whereas Leptopilinaboulardi is a specialist. Recent studies have revealed that these wasps have developed different strategies to defeat the immune defense of Drosophila. Whole-genome sequencing (WGS) analysis has been performed, and distinct sets of venom proteins have been reported.

In this study, we focussed on Asobara japonica Belokobylskij (Hymenoptera: Braconidae). Asobarajaponica is a typical generalist parasitoid wasp that lays eggs in the larvae of various Drosophila species., Similar to the genus Leptopilina, an adult female A. japonica wasp oviposits an egg into the host larva, enabling a wasp larva to hatch inside the host. The parasitoid wasp larva grows with the host fly larva until pupariation. Inside the host pupa, the wasp kills the host after pupariation and eventually ecloses from the host pupal case. It takes 15–17 days to become an adult wasp after oviposition. Notably, A. japonica has a high parasitic success rate (>90%), implying that the host immune defense system hardly affects wasp development. This feature enables us to consistently study the developmental process of parasitoid wasps, in contrast to studying the ‘tug-of-war’ relationship between host immune defense response and parasitic strategy. Moreover, it has been reported that the venom of A. japonica has a severe deleterious effect on the host larvae,: Interruption of their oviposition behaviour before egg laying causes high larval mortality in host Drosophila species. This suggests that the presence of toxic venom components as well as their neutralizer components underlie the unique parasitic strategy in A. japonica, which is distinct from the venom components of other parasitoid wasps.

Based on host specificity and the death of host animals, many parasitic species are used as biological control agents for the management of pest species. Recent studies have evaluated the parasitism of A. japonica, because it parasitizes Drosophila suzukii (Matsumura), commonly known as spotted wing Drosophila, which is a serious economic threat to the production of fruits worldwide., Considering the deleterious effect of A. japonica on the hosts,, it is crucial to understand the molecular mechanisms of parasitism in order to develop novel insecticide seed strategies.

Here, we performed WGS analysis of A. japonica and devised a protocol for RNA interference (RNAi) with A. japonica. To obtain a homogenous genome sequence, we generated the clonal strain ‘Genome #3’ from a single female wasp and successfully constructed the highly qualified WGS platform for gene annotation, compared with the previous study. The well-organized WGS platform allows us to conduct genetic approaches from both hosts and parasites in a bidirectional manner.

2. Materials and methods

2.1 Biological materials

A parthenogenetic strain of the endoparasitoid wasp A.japonica ‘Tokyo (TK)’ was originally collected near Tokyo by Dr. Masato Kimura and kindly provided by Dr. Kazuo Takahashi at Okayama University. The clonal strain ‘Genome #3’ was established from a single female in the Tokyo strain. In D.melanogaster, the Oregon R (OR) strain was used as the wild-type and host in this study. The e mutant (DGRC #106436) was obtained from KYOTO Stock Center (DGRC) at Kyoto Institute of Technology. A.japonica and D. melanogaster were reared on standard agar-cornmeal medium at 25°C under a 12:12-h light/dark cycle.

2.2 Genome preparation

In the clonal strain ‘Genome #3’, 200 Wolbachia-free wasps were collected from host fly larvae reared on a standard cornmeal diet supplemented with tetracycline (4 mg/2 g food). Wasps were flash-frozen in liquid nitrogen and stored at −80°C. Genomic DNA was extracted using the genomic DNA buffer set and Genomic-tip 100/G (QIAGEN). After the genomic DNA samples were resuspended in 1× Tris-EDTA buffer (pH 8.0), the concentration of DNA was measured using a Nanodrop and Qubit dsDNA BR Assay kit (Thermo Fisher Scientific). The quality of the samples was examined by DNA electrophoresis. A total of 38 µg of 200 kb DNA was used for WGS.

2.3 Genome sequencing and assembly

Genome sequencing was performed using Illumina HiSeq 2500, for paired-end short reads, and PacBio Sequel II platforms. Genome size and heterozygosity were estimated using the GenomeScope 2.0., The primary assembly was performed using FALCON-UNZIP with default settings, and the assembled contigs were polished using Arrow and Pilon. Completeness of the assembly was assessed using BUSCO v5.2.2 with the ‘insecta_odb10’ and ‘hymenoptera_odb10’ datasets. For comparison, we performed BUSCO analysis with other Hymenopteran species: Apis mellifera, Orussus abietinus, Trichogramma pretiosum, Copidosoma floridanum, Nasonia vitripennis, Ceratosolen solmsi, Diachasma alloeum, Fopius arisanus, Microplitis demolitor, Macrocentrus cingulum, Cotesia typhae, Cotesia glomerata, Chelonus insularis, and Aphidius gifuensis. The references of the genomes of all species are listed in Supplementary Fig. S1.

2.4 RNA sequencing

Total RNA was extracted from whole bodies of 200 female wasps, venom glands of 300 female wasps at 1 day after eclosion, and venom glands of 300 female wasps at 10 days after eclosion. Samples were dissected and ground with pestles in RNAiso Plus reagent (TaKaRa). After DNase treatment, the concentration of RNA was measured using a Nanodrop and Qubit double-stranded RNA (dsRNA) BR Assay kit (Thermo Fisher Scientific). Sequencing libraries were prepared by means of polyA selection step with TruSeq Stranded mRNA Library Prep kit (Illumina). Then, libraries were sequenced on an Illumina HiSeq 2500 (101 cycles) for paired-end short reads with a length of 100 bp. After filtering the low-quality and adapter sequences, 82,282,218, 84,554,634, and 85,120,902 RNA-Seq reads were obtained, respectively.

2.5 Gene annotation

Protein-coding gene prediction was conducted using an in-house annotation pipeline based on a combination of RNA-seq-based prediction, homology-based prediction, and ab initio prediction methods. To map the first RNA-seq-based prediction, RNA-seq reads were mapped to assembled genomes using HiSat2 (version 2.1.0), and gene structure was predicted using StringTie (version 1.3.4d). For de novo first RNA-seq-based prediction, Trinity and Oases were used. Assembled RNA contigs were spliced and mapped by GMAP and the open reading frame (ORF) was predicted using ORFfinder. For homology-based predictions, protein sequences from four species (D. alloeum, F. arisanus, M. demolitor, and M. cingulum) were aligned to the assembled genomes using Spaln. For ab initio gene prediction, Augustus and SNAP were trained with RNA-seq-based predicted genes and then used. Moreover, functional annotation was performed using BLASTp for homology signatures against NCBI databases and InterProScan for Pfam, and the gene ontology terms associated with the proteins were retrieved using the InterPro ID. The subcellular localization was predicted using SignalP.,

2.6 Phylogenetic analysis

To compare orthologues, we used protein sequences of 14 hymenopteran species available from three repositories. Protein sequences of O. abietinus, D. alloeum, F. arisanus, T. pretiosum, M. cingulum, M. demolitor, C. floridanum, N. vitripennis, and C. solmsi were downloaded from WaspBase (http://www.insect-genome.com/waspbase/); A. mellifera, C. typhae, C. insularis, and A. gifuensis from NCBI (https://www.ncbi.nlm.nih.gov); and C. glomerata from Figshare (https://figshare.com). Subsequently, the identification of orthologous groups of proteins was performed. First, protein sequences from all species were grouped into gene families and 2,381 orthologue groups with a one-to-one relationship across all species were extracted using Proteinortho. For each group, multiple alignments were performed with MAFFT and sites containing gaps (‘–’) or ambiguous characters (‘X’) were excluded. All alignments were concatenated, and 1,110,464 amino acid sites were used for the phylogenetic analysis. A phylogenetic tree was constructed using RAxML (version 8.2.12). Here, we applied the JTT substitution matrix with a gamma model of rate heterogeneity (-m PROTGAMMAJTT).

2.7 Data availability

The raw sequence data generated on the Illumina and PacBio platforms were deposited in the DNA Databank of Japan (DDBJ, BioProject ID: PRJDB12783). The accession numbers of genome assembled scaffold data are BQMI01000001–BQMI01000601. Gene annotation data are available upon request.

2.8 Observation of the developmental progression of A. japonica

Wild-type fly eggs were laid for 2 h on a grape plate with yeast paste to synchronize host development. After 24 h, the first instar (L1) larvae were transferred into vials containing standard cornmeal food. At 96 h after egg-laying, the third instar (L3) larvae and adult wasps were transferred into the bulk infection arena (Fig. 3A). Wasps were allowed to freely lay eggs in L3 host larvae for 1 h and were then removed from the arena to synchronize the timing of oviposition. At a certain time point after infection, a wasp larva was carefully removed from the host pupal case with forceps under a dissection microscope and placed on a 2% agar plate. When the wasp larvae were removed from the host pupal case at 4 days post-infection (dpi), development proceeded on the agar plate and wasp adults eclosed normally.

Figure 3

Schematic representation of the RNAi protocol during A. japonica development. (A) A standard setup of a bulk infection arena for wasp parasitism. (B) Representative images of developing wasps on an agar plate. dpi, days post-infection. Scale bar, 1 mm.

2.9 Identification of the A. japonica ebony gene and putative venom genes

We identified the A. japonica ebony gene (gene003270, ebony-Aj) as an orthologue of the D. melanogaster ebony (ebony-Dm) gene. Ebony encodes the enzyme N-β-alanyl dopamine (NBAD) synthetase, which converts dopamine to NBAD. Loss of this enzymatic function increases dark melanin and decreases tan pigment, resulting in a dark coloured body phenotype. Thus, we chose ebony-Aj to be the target of RNAi for the visible loss-of-function phenotype in the first trial. The D. melanogaster amino acid sequence of the ebony gene was obtained from FlyBase (http://flybase.org/). The A. japonica orthologue of ebony (ebony-Aj) was identified by tBLASTn analysis using the amino acid sequence of Ebony-Dm as a query.

To identify putative venom genes, the fragments per kilobase of transcript per million mapped reads (FPKM) values were calculated with RNA-seq data of the venom gland sample (1 day after eclosion). Adapter sequences were trimmed from raw RNA-seq reads using Trimmomatic (version 0.39). RNA-seq reads were mapped to assembled genomes using HiSat2 (version 2.1.0). Mapped read data were sorted, merged and counted using SAMtools (version 1.9) and StringTie (version 2.0.6). FPKM values were calculated with R (version 3.6.1) and Ballgown (version 2.18.0).

2.10 Quantitative reverse transcription PCR

To confirm if these putative venom genes are strongly expressed in the venom gland, we performed quantitative reverse transcription PCR (qRT-PCR) analysis of the top five genes with high FPKM values (Table 7). The venom glands and the other carcass from 20 wasp adults were dissected and total RNA was extracted using RNAiso Plus reagent (TaKaRa). cDNA was prepared with ReverTra Ace qPCR RT Master Mix with gDNA Remover (TOYOBO). Quantitative PCR was conducted using the Universal SYBR Select Master Mix (Applied Biosystems) with a Thermal Cycler Dice TP800 system (TaKaRa). For the relative quantification of expression, the A. japonica orthologue of Ribosomal protein L32 (RpL32-Aj) was identified by tBLASTn analysis using amino acid sequence of D. melanogaster RpL32. The target gene expression level was normalized to RpL32 and then relative fold changes were calculated by the delta-delta Ct method. The mean values were from three independent experiments. The primers used for qRT-PCR are described in Table 1.

Table 7

Top five genes with high FPKM values in the venom gland sample

	Gene ID	FPKM value	Relative expression level (venom gland/other carcase)
			Mean value (n = 3)
1	gene006533	403151.38	110,082.6
2	gene010987	395514.72	N.A.
3	gene010975	113373.65	1,401,477.8
4	gene003054	112626.62	1,185,641.8
5	gene000789	34218.91	258,595.8

N.A. (not available): The coding sequence of gene010987 was an exact match with that of gene006533 and the gene locus of gene010987 was just next to gene006533. qRT-PCR was not performed for gene010987.

Table 1

Primer sequences used in this study

Primer name	Sequence (5′> 3′)	Purpose
Ajebony_ex8_Fwd	TATTTTAACCGAAAGTTTCTACTTGAAAGCTG	dsRNA template
Ajebony_ex8_Rev	CTGTGTAAGAGGACTCGTATTGGTG	dsRNA template
AjRpL32_qPCR_Fwd	CCCGTCACATGCTTCCTACT	qRT-PCR
AjRpL32_qPCR_Rev	GAATTTGCGATTCTGCATCA	qRT-PCR
Aj_gene003054_ex1_Fwd	ATGTACTCCCACTGTAGGTTCCAAG	CDS cloning
Aj_gene003054_ex3_Rev	TTACTTCCTCCCCGTAAGCGCTC	CDS cloning, dsRNA template
Aj_gene003054_RNAi_Fwd	CAAGAAAGCAAACGGCAAACCTTGG	dsRNA template
T7_Aj_gene003054_RNAi_Fwd	GGATCCTAATACGACTCACTATAGCAAGAAAGCAAACGGCAAACCTTGG	dsRNA template
T7_Aj_gene003054_ex3_Rev	GGATCCTAATACGACTCACTATAGGTTACTTCCTCCCCGTAAGCGCTC	dsRNA template
qPCR_gene006533_Fwd	CATCGGAACTACAGGGCATT	qRT-PCR
qPCR_gene006533_Rev	TGCCAATGTCTTCACACTCC	qRT-PCR
qPCR_gene003054_Fwd	CTGATTGTCGTGCTCGGTTA	qRT-PCR
qPCR_gene003054_Rev	CCACCCTGAGGATGTGTTTC	qRT-PCR
qPCR_gene010975_Fwd	GTATCTTCGGGATGCTCTGC	qRT-PCR
qPCR_gene010975_Rev	CCCTCCGCTAACTCACACAT	qRT-PCR
qPCR_gene010975_Fwd2	CCTGATAATCGTCGGTATCTTC	qRT-PCR
qPCR_gene010975_Rev2	TGCCCACTGTTCCTGACATC	qRT-PCR
qPCR_gene000789_Fwd	TTCTGAGACAGAGCCCGAAT	qRT-PCR
qPCR_gene000789_Rev	GGAATATGCAGTGGGTCGTC	qRT-PCR

Underline indicates the T7 promoter sequence. dsRNA, double-strand RNA; qRT-PCR, quantitative reverse transcription PCR; CDS, coding sequence.

2.11 Primer design for RNAi target regions

To conduct ebony-Aj RNAi experiments, we chose the exon 8 sequence as an RNAi target region because it is the longest single exon among the ebony-Aj (gene003270) gene region. For the putative venom gene gene003054, the 151 bp region of the coding sequence was chosen as the RNAi target region. RNAi target regions were distinct from the qPCR target regions. For qPCR, each primer was designed using Primer3Plus. The primer sequences are listed in Table 1.

2.12 dsRNA synthesis and injection

For dsRNA synthesis, a gene-specific sequence was used as a template for transcription. Template cDNAs were obtained from total RNA extracted from the whole body of A. japonica using RNAiso Plus reagent (TaKaRa) and PrimeScript™ Reverse Transcriptase (TaKaRa). PCR was performed using KOD Plus Neo (TOYOBO) and the amplified specific DNA fragments were inserted into the SmaI site of the pBluescript SK (-) plasmid with Ligation high Ver. 2 (TOYOBO). These plasmids were digested with NotI and EcoRI restriction enzymes, and specific DNA fragments were ligated into NotI and EcoRI cut pBluescript KS (+) in the opposite direction to the T7 promoter sequence. Finally, these two complementary DNA fragments were digested with NotI and EcoRI to generate linearized DNA templates in the T7 RiboMAX™ Express RNAi system (Promega) following the manufacturer’s instructions. The Green fluorescent protein (GFP)-coding sequence containing the p-GEM easy plasmid was a gift from Dr. Manabu Kamimura (National Agriculture and Food Research Organization, NARO) and was used for the synthesis of a control GFP dsRNA.

dsRNA microinjection was performed with a glass needle made from a glass capillary with filament (GD-1, NARISHIGE) using a puller (PC-10, NARISHIGE) and polished with a micro grinder (EG-401, NARISHIGE). The dsRNA solution was injected into the abdomen of A. japonica larvae, after the host pupal case was carefully removed and then placed on a 2% agar plate. For ebony-Aj RNAi, dsRNA was concentrated to 10 μg/μl by ethanol precipitation for microinjection. In the case of the putative venom gene-RNAi, dsRNA solution was diluted to 1,000, 500, and 250 ng/μl with ultrapure water. Although we could not control the exact amount of synthesized dsRNA for microinjection because of limitations imposed by our injection apparatus, we could estimate the required amount of 100–200 nl.

2.13 Evaluation of body colour

The dorsal side of the wasp bodies was imaged using a digital camera (Leica MC120 HD) attached to a dissection microscope under constant lighting conditions. A square area of the wasp thorax was extracted to quantify the body colour of individual wasps. To minimize the effects of highlights and shadows in the images, we selected an area devoid of pixels with overexposure brightness and underexposure darkness. Then, 8-bit integer RGB values (0 ≤ R, G, B ≤ 255) were measured and averaged using Fiji software (National Institutes of Health), and the redness index was calculated by Σ(R—mean[R, G, B])/n (n = the number of pixels).

3. Results and discussion

3.1 Generation of a clonal strain derived from a single female wasp

To obtain a homogenous genome sequence with high quality, it is essential to collect enough amount of uniform genome DNA from individual samples. In A. japonica, the Wolbachia-infected strain ‘Tokyo (TK)’ is available, which produces only female wasps from unfertilized eggs by Wolbachia-induced parthenogenesis. Therefore, we generated a clonal strain ‘Genome #3’ from a single female wasp in TK strain (Fig. 1). After Genome #3 was established, Wolbachia was removed from these wasp clones to avoid contamination with Wolbachia DNA of the wasp genomic DNA samples. When host fly larvae were reared on a standard diet supplemented with tetracycline, parasitoid wasp eggs and larvae were exposed to tetracycline through the haemolymph of the hosts and Wolbachia was eliminated from wasp bodies., Under these conditions, eclosed female wasps produced not only female but also male offspring in the next generation. We confirmed the removal of Wolbachia at 4 mg tetracycline in 2 g fly food, and collected 200 female wasps with an identical genetic background for genomic DNA extraction (Supplementary Fig. S2A).

Figure 1

Generation of a clonal strain ‘Genome #3’ from a single female wasp. To obtain the high-quality homogenous genome sequence, we generated a clonal strain from a single female wasp derived from the Wolbachia-infected parthenogenetic strain ‘Tokyo’. After the clonal strain ‘Genome #3’ was established, we removed Wolbachia from these wasp clones by rearing host fly larvae on the standard food supplemented with tetracycline. 200 Wolbachia-free wasps were collected for genome DNA extraction. At the same time, we confirmed the removal of Wolbachia by the appearance of male wasp offspring in the next generation.

3.2 Genome assembly and assessment of genome completeness

We sequenced DNA extracted from adult female A. japonica wasps with an identical genetic background. PacBio Sequel of 27.7 Gb on 1,627,155 reads (average read length: 16 kb) and Illumina Hiseq2500 reads of 64.3 Gb on 257,165,308 read-pairs (read length: 250 bp, average insert size: 631 bp) were obtained (Table 2). From the short reads, the genome size was estimated to be 322 Mb., The heterozygosity was estimated to be 0.132%, indicating that a homogenous genome sequence was successfully obtained (Supplementary Fig. S2B). This value reflects the clonal population of Genome #3. Assembly using FALCON-UNZIP resulted in 601 primary contigs with a total length of 322.2 Mb and an N50 of 2.64 Mb (Table 2). This assembled size is highly consistent with the estimated genome size and comparable to that of other strains of A. japonica, which was recently reported to be 281–304 Mb (Kagoshima, asexual strain) and 280–284 Mb (Amami-Oshima, sexual strain). The total length of our WGS was 322.2 Mb, which is longer than that of previous study (270.9 Mb), indicating that our WGS covers broader, unspecified genomic regions than the previous study (Table 3). Moreover, the genome size is comparable to that of other Hymenopteran species (Table 3).

Table 2

Determination of the WGS of A. japonica

(1) PacBio Sequel
Total number of bases (bp)	Total reads	Max read length (bp)	N50 (bp)	Average read length (bp)
27,657,055,104	1,627,155	213,716	29,875	16,997

N50 is the shortest contig length that needs to be included for covering 50% of the genome.

Assembly: Falcon (v0.7) + FalconUnzip (v0.4.0) error correction: Arrow (SMRT Link v6.0.0) + Pilon (v1.22).

Table 3

Genome assembly values among Hymenopteran species

(1) Comparison of statistics of genome
	Total length (Mb)	Number of scaffolds	Average of scaffold (kb)	Longest scaffold (Mb)	N50 (Mb)	L50	GAP (%)
A. japonica (this study)	322.2	601	536.1	7.9	2.6	38	0
A. japonica 18	270.9	774	350	5.8	1.7	47	0
A. mellifera a	225.3	177	1,272.6	27.8	13.6	7	0.6
O. abietinus	201.2	936	215	6.44	2.4	28	7.3
T. pretiosum	195.1	357	546.5	8.8	3.7	19	7.7
C. floridanum	555.0	5,445	101.9	8.1	1.0	153	18.0
N. vitripennis	295.8	6,098	48.5	33.6	0.9	21	19.3
C. solmsi	277.1	2,457	112.8	27.4	9.6	10	0.5

The Bold values were discriminated from the other values that have been already reported in other references.

Refseq is used for analyses. Refseq is a curated genome data by NCBI using genome data which is uploaded by users.

The completeness of the genomic assembly data was evaluated by detecting a set of unique single-copy genes in the primary sequence with Benchmarking Universal Single-Copy Orthologues (BUSCO, v5.2.2). Of the 1,353 BUSCOs in the insect dataset ‘insecta_odb10’, 1,337 (98.8%) complete BUSCOs were detected in the assembly (Table 3). The number of complete and single-copy BUSCOs was 1,326 (98.0%). In addition, we also performed BUSCO analysis using the Hymenoptera (‘hymenoptera_odb10’; 5,991 BUSCOs) dataset and identified 93.8% complete BUSCOs in our assembly (Table 3). This rate is sufficiently high and comparable to that of the previous study (92.5%) and other species, including A. mellifera (97.7%) and N. vitripennis (89.7%). In the family of Braconidae, complete BUSCOs were calculated in D. alloeum (95.3%), F. arisanus (94.5%), M. demolitor (94.4%), M. cingulum (93.5%), C. typhae (93.2%), C. glomerata (93.6%), C. insularis (93.5%), and A. gifuensis (91.3%) (Table 4). Taken together, these results indicate that the generated assembly was of high quality for gene annotation in A. japonica.

Table 4

Genome assembly values among Blaconidae species

(1) Statistics of genome
	Total length (Mb)	Number of scaffolds	Average of scaffold (kb)	Longest scaffold (Mb)	N50 (Mb)	L50	GAP content (%)
A. japonica (this study)	322.2	601	536.1	7.9	2.6	38	0
A. japonica 18	270.9	774	350	5.8	1.7	47	0
D. alloeum	388.8	3,968	98.0	6.6	0.6	128	6.2
F. arisanus	153.6	1,042	147.4	5.5	1.0	49	8.2
M. demolitor	241.2	1,794	134.4	7.1	1.1	50	14.7
M. cingulum	132.4	5,696	23.2	1.4	0.2	179	3.0
C. typhae	186.7	72	2,592.5	20.0	6.8	8.0	0
C. glomerata	288.8	3,354	86.1	37.5	27.8	5.0	0
C. insularis a	135.7	455	298.3	4.7	1.2	33.0	0
A. gifuensis	156.9	24	6,539.4	32.1	27.5	3.0	0

Refseq is used for analyses. Refseq is a curated genome data by NCBI using genome data which is uploaded by users.

3.3 Gene prediction, annotation, and phylogenetic tree analysis with hymenopteran species

Protein-coding gene prediction was conducted using an in-house annotation pipeline based on a combination of RNA-seq-based prediction, homology-based prediction, and ab initio prediction methods. The final gene model set consisted of 12,508 genes with a mean coding sequence (CDS) length of 1,594 bp (Table 5). The statistics of gene annotation data were comparable to other Hymenopteran species. Moreover, the completeness of the gene annotation data was evaluated by BUSCO in the insect dataset and the Hymenoptera dataset: 12,283 (98.2%) and 11,933 (95.4%) complete BUSCOs were detected in the annotated genes, respectively. We also compared these rates with other Braconidae species (Table 6). Overall, these results strongly support that our gene annotation data is reliable for the identification of genes in A. japonica.

Table 5

Comparison of gene annotations between A. japonica and other Hymenopteran species

(1) Statistics of gene annotation
	Number of genes	Number of exons per gene	Number of single exon gene	Total exon length (Mbp)	Mean exon length (bp)	Mean CDS length (bp)	Total intron length (Mbp)	Mean intron length (bp)	GT-AG splicingsite (%)
A. japonica	12,508	5.7	608	19.9	282.2	1,594.1	103.2	1,774.5	99.4
A. mellifera a	9,922	7.1	503	18.1	257.4	1,824.1	102.6	1,699.7	98.6
O. abietinus	9,996	7.0	533	18.1	260.0	1,810.2	77.0	1,290.9	98.6
T. pretiosum	12,422	5.8	1,238	21.0	290.5	1,689.5	67.9	1,134.6	98.3
C. floridanum	11,907	6.2	989	19.2	259.2	1,612.8	195.7	3,147.1	98.0
N. vitripennis	13,185	6.1	1,226	21.0	259.2	1,590.6	86.9	1,283.4	98.0
C. solmsi	9,823	6.9	601	16.5	244.9	1,683.8	106.0	1,837.2	96.5

Refseq is used for analyses. Refseq is a curated genome data by NCBI using genome data which is uploaded by users.

Table 6

Comparison of gene annotations between A. japonica and other Blaconidae species

(1) Statistics of Gene annotation
	Number of genes	Number of exons per gene	Number of single exon gene	Total exon length (Mbp)	Mean exon length (bp)	Mean CDS length (bp)	Total intron length (Mbp)	Mean intron length (bp)	GT-AG splicing site (%)
A. japonica	12,508	5.7	608	19.9	282.2	1,594.1	103.2	1774.5	99.4
D. alloeum	12,837	4.8	1,759	21.1	338.7	1,640.8	124.4	2,520.0	98.8
F. arisanus	10,991	5.0	819	18.5	334.9	1,686.1	53.4	1,204.5	98.4
M. demolitor	12,144	5.3	1,339	20.4	318.9	1,676.5	79.1	1,529.7	98.2
M. cingulum	11,993	5.1	1,218	18.5	300.5	1,541.5	19.0	384.1	100.0
C. typhae	8,591	5.0	483	12.8	299.9	1,487.4	20.4	598.6	99.4
C. glomerata	19,218	4.2	820	22.8	284.4	1,187.8	45.6	746.7	99.7
C. insularis a	10,548	5.6	826	18.7	315.7	1,773.4	45.2	928.6	98.8
A. gifuensis	11,504	5.4	60	19.5	316.3	1,695.3	48.6	969.6	99.3

Refseq is used for analyses. Refseq is a curated genome data by NCBI using genome data which is uploaded by users.

Single-copy orthologous genes were determined from 10 species (A. japonica, D. alloeum, F. arisanus, M. demolitor, M. cinglum, C. typhae, C. glomerata, C. insularis, A. gifuensis, and A. mellifera) and 3,263 groups were extracted with a one-to-one relationship across all species. For each group, multiple alignments were performed using the MAFFT. A phylogenetic tree was constructed using RAxML (version 8.2.12). The orthologues of A. mellifera were used as the outer group. A total of 1,000 replicates were used for bootstrap analysis (Fig. 2). Our results are consistent with the phylogeny based on mitochondrial DNA, as previously reported.,

Figure 2

Phylogenetic analysis of A. japonica and other Hymenopteran species. (A) The phylogenetic tree was generated with 2,381 pairs of single-copy orthologues between A. japonica and 14 related Hymenoptera species. (B) The phylogenetic tree was generated with 3,263 pairs of single-copy orthologues between A. japonica and 8 related species in the family Braconidae. Apis mellifera was used as an out group. All nodes have 100% bootstrap support after 1,000 replications.

3.4 RNAi technique protocol with A. japonica

To the best of our knowledge, there have been no previous studies on genetic manipulation in A. japonica. To perform functional analyses of A. japonica genes, we devised a protocol for dsRNA microinjection-based RNAi. This procedure is a post-transcriptional gene silencing mechanism that applies exogenous dsRNA to wasp bodies. dsRNAs are cleaved into small interfering RNA fragments that bind to complementary endogenous mRNAs. Target mRNAs result in degradation such that the corresponding gene products are decreased. The effects of RNAi have already been reported in various Hymenoptera species, including the ectoparasitoid wasp N.vitripennis,, the endoparasitoid wasp Microplitis,, and L. boulardi. Therefore, we reasoned that these methods could be applied to experiments with A. japonica.

As in L. boulardi, A. japonica larvae undergo metamorphosis inside the host puparium, which makes it difficult to conduct dsRNA microinjection experiments. To examine whether the development of A. japonica progresses outside the host puparium, we set up a bulk infection arena to synchronize the timing of oviposition (Fig. 3A). In addition, third-instar (L3) larvae (at 96 h after egg-laying) were used as hosts. At a certain time point after infection, we carefully removed the host pupal case with forceps and observed wasp development on the 2% agar plate (Fig. 3B). When the wasp larvae were removed from the host pupal case at 4 dpi, development proceeded on the agar plate and wasp adults eclosed normally. During this time, the wasp larva completely consumed the fly host body. We recommend that the appropriate timing for dsRNA injection is after 4 dpi.

For the RNAi experiment, we chose the ebony gene to be the target of RNAi for the first trial. Ebony encodes the enzyme NBAD synthetase, which converts dopamine to NBAD. Loss of this enzymatic function increases dark melanin and decreases tan pigment, resulting in a dark body colour phenotype. We identified gene003270 (ebony-Aj) as an orthologue of the D. melanogaster ebony (ebony-Dm) gene from the annotated A. japonica gene set using BLAST.

It appeared that eye and body colour pigmentation started from 7 to 9 dpi (Fig. 3B). We removed the pupal case and injected dsRNA into A. japonica at 6 dpi, prior to body colour pigmentation. Therefore, the body colour of ebony-Aj dsRNA-injected wasps became darker than that of the control GFP dsRNA-injected wasp (Fig.4A and B and Supplementary Fig. S3). To quantitatively evaluate the body colour, we calculated the redness index from the captured images (Fig. 4C). In our analysis, the redness index was significantly decreased in ebony-Aj dsRNA-injected wasps, which is similar to the body colour difference between Drosophila wild-type and ebony mutant animals (Supplementary Fig. S4). From this result, we concluded that the function of Ebony was successfully silenced in A. japonica using RNAi.

Figure 4

RNAi effects of RNAi wasps were evaluated by phenotype analysis and qPCR. (A) Representative images of wasp adults being injected with GFP dsRNA as a control. A square indicates the region of interest (ROI) for RGB value measurements. (B) Representative images of wasp adults being injected with ebony-Aj dsRNA as ebony RNAi. (C) Quantitative evaluation of wasp body colour. The Y-axis indicates the redness of the ROI. ***P < 0.005 from Student’s t-test. n = 18 (GFP dsRNA), 21 (ebony-Aj dsRNA). (D, E) The relative expression levels of gene003054 (D) and gene010975 (E) in GFP-RNAi or gene003054-RNAi wasps were quantified using the delta-delta Ct method. The expression of RpL32-Aj was used to normalize the values. All values represent the means ± SD with all data points (n = 3). *P < 0.05, **P < 0.01, and ***P < 0.005 from Student’s t-test. n.s., non-significant (P > 0.05).

3.5 Knock down of the putative venom gene in A. japonica

After confirming the RNAi effect by ebony-Aj-dsRNA microinjection, we applied this method to knock down the putative venom gene in A. japonica. Our RNA-seq analysis indicates that a group of genes are strongly expressed in the venom gland, which is likely to be involved in venom production. We listed these putative venom genes in descending order of FPKM values and examined the expression of the top five genes (gene006533; gene010987; gene010975; gene003054; gene000789) by qRT-PCR (Table 7). Consistent with RNA-seq data, these genes were exclusively expressed in the venom gland compared with other tissues. These gene products showed little homology signatures against NCBI database and few functional domains against Pfam and SMART, suggesting that they are specific to A. japonica. Among them, gene003054 encodes a putative transmembrane protein, and we chose it as a target of RNAi.

dsRNA microinjection was performed with wasp larvae at 6 dpi. GFP-dsRNA was used as the negative control. To determine the appropriate concentration of dsRNA, we used 1,000, 500, and 250 ng/μl dsRNA-containing ultrapure water for microinjection. When dsRNA-injected wasps were eclosed, total RNA was extracted from these wasps, and the relative expression level of each gene was measured by qRT-PCR. In sharp contrast to GFP-dsRNA injection, the expression levels of the target gene003054 were drastically reduced in all three concentrations of dsRNA-injected wasps (Fig. 4D). To test the specificity of RNAi, we also examined the expression of gene010975 (Fig. 4E). The expression did not decrease in 1,000, 500, and 250 ng/µl dsRNA injection in gene003054–RNAi wasps. These results suggest that dsRNA injection-based RNAi specifically suppressed the expression of the target gene (Fig. 4D and E). In addition, 250 ng/μl of dsRNA was sufficient for effective RNAi against putative venom genes. After eclosion, the dsRNA-injected wasps were viable but less tolerant to the behavioural assay. It is possible that microinjection or rearing conditions may affect wasps with a weak constitution. Taken together, we confirmed that the RNAi method provides an efficient way to silence putative venom genes in A. japonica.

4. Conclusions

Here, we performed WGS analysis of A. japonica, annotating a set of 12,508 genes. These genes include highly conserved sets of genes in Hymenoptera, indicating the high reliability of genome data information. The qualified WGS platform provides us with the references needed to identify venom components for parasitizing Drosophilidae and the diversity of species during coevolution. Considering the virulence effects of venoms and the specificity of host insects, including D. suzukii, venom genes and their active products are crucial for developing insecticides in agricultural sciences. Moreover, we succeeded in devising a protocol for RNAi techniques using A. japonica. Given the genetic tools available in the host fly D. melanogaster, the development of genetic tools in the parasitoid wasp A. japonica enables a bidirectional approach to examine host-parasite interactions. Based on the WGS and functional analyses, we confirmed that it is possible to explore these molecular mechanisms in non-model organisms.

Supplementary data

Supplementary data are available at DNARES online.

Author contributions

T.Ka., M.K., H.M., T.Ku., A.K., and Y.S.-N. conducted the experiments. A.T., H.T., T.I., and T.Ka. performed WGS and RNA-seq analyses. T.Ka., Y.S.-N, H.T., and R.N. wrote the article with feedback from all the co-authors.

Click here for additional data file.