In northern China, in addition to L. dispar, B. lasus is also an important pupal parasitoid of H. cunea, for which the parasitism ratio is approximately 1.06–3.39% in the field (Yang et al., 2001). To acquire B. lasus adults, we collected the pupae of H. cunea, which may be parasitized by B. lasus and other parasitoid species (e.g., Coccygomimus disparis Viereck; Chouioia cunea Yang) from a field in Xuzhou City, Jiangsu Province, China, in March 2016. After collection, we isolated the pupae individually in polyethylene tubes (height: 7.5 cm; diameter: 1 cm) whose openings were covered by a cotton ball and incubated them at a temperature of 28 ± 0.5°C, a relative humidity (RH) of 70 ± 5% and a photoperiod of 14 L:10 D. We observed and selected B. lasus after adult eclosion.

For the transcriptomic experiment, only 1-day-old B. lasus adults were selected, and the sex was determined under a microscope (Leica M205A, Germany). Then, five adults of the same sex were pooled into a plastic tube (1.5 ml), snap-frozen in liquid nitrogen, and transferred to a –80°C freezer for long-term storage. RNA from each sample group (whole bodies of female and male adults) was extracted with TRIzol reagent (Invitrogen; United States). Each group had three replicates. The quality of the isolated RNA was assessed using a NanoDrop (Thermo Fisher Scientific NanoDrop 2000, United States), and the A260/280 values were all above 2.0. A total of 3 μg total RNA from each sample was converted into cDNA using the NEBNext^® Ultra^TM RNA Library Prep Kit for Illumina^® (NEB, United States). In total, six cDNA libraries were constructed and subsequently sequenced with the Illumina HiSeq 2000 platform by Beijing Biomarker Technologies Co., Ltd, resulting in raw reads. Raw sequence data generated were deposited into Sequence Read Archive (SRA) database of NCBI with the accession no. PRJNA513855. Clean reads were obtained by removing reads containing adapter, poly-N reads and low-quality reads from the raw data using FASTX-Toolkit¹, and these clean reads were used for further analysis. Then, transcriptome assembly was performed using Trinity (v2.5.1) with the default parameters (Grabherr et al., 2011). For functional annotation, pooled assembled unigenes were searched using BLASTX (v2.2.31) against five public databases, Clusters of Orthologous Groups (COG), Swiss-Prot, NCBI non-redundant protein sequences (nr), KEGG Ortholog database (KO) and GO, with an E-value cutoff of 10^-5. Using our assembled transcriptome as a reference, we identified putative genes expressed in males and females by RSEM (Li and Dewey, 2011), using the reads per kb per million reads (RPKM) method. Genes with at least 2-fold changes (i.e., log₂∣FC∣≥ 1) and a false discovery rate [FDR] < 0.01 as found by DESeq R package (1.10.1) were considered differentially expressed. The GOseq R package (Young et al., 2010) and KOBAS software (Mao et al., 2005) were used to implement the statistical enrichment of differentially expressed genes (DEGs) in the GO and KEGG pathways, respectively, and an adjusted Q-value <0.05 was chosen as the significance cutoff.

Based on transcriptomic data, a gene of transient receptor potential (trp) involved in the phototransduction pathway enriched only in males (ko: 04745; Supplementary Figure S1-d), trp (Leung et al., 2000), was down-regulated in females, which may lead to a reduction in light response (Leung et al., 2000; Popescu et al., 2006). Therefore, we checked this result at the mRNA expression and behavioral levels.

Prior to analysis, the raw data were tested for normality and homogeneity of variances with the Kolmogorov-Smirnov test and Levene’s test, respectively, and the data were log-transformed if necessary. The qRT-PCR data comparing gene expression in females and males were analyzed with the independent t-test. In phototaxis assays, the preferences for light and dark were analyzed using sign tests, and the differences in female and male phototaxis were analyzed by the chi-square test. All analyses were performed using SPSS v.20 (IBM SPSS, Armonk, NY, United States).

All high-quality reads (101,945,678) from the six samples were pooled and assembled by using Trinity with the default parameters, and a total of 254,656 transcripts with lengths longer than 200 bp were generated. The N50 size was 2706 bp with 57,605 sequences longer than 1 kb. We chose the longest isoform of each gene to construct the unigene set. After isoforms were considered, these assembled transcripts were predicted to be produced from a total of 164,709 unigenes. The N50 size of the unigenes was approximately 814 bp, and their mean length was 572.08 bp (Supplementary Table S2). For annotation, the pooled assembled unigenes were searched using blastx against five public databases with an E-value cutoff of 10^-5. A total of 37,453 unigenes were successfully annotated, as shown in Table 1, including 17,248 genes in GO, 13,491 genes in COG, 35,427 genes in nr, 18,195 genes in Swiss-Prot, and 15,133 genes in KEGG.

In the GO analysis, 17,248 unigenes were successfully annotated and classified into three major GO categories: molecular function (MF), cell component (CC), and biological processes (BP), then assigned to 56 subcategories based on GO level 2. The dominant subcategories for the classified genes were catalytic activity and binding for the MF category; cell and cell part for the CC category; and metabolic process, cellular process, and single-organism process for the BP category (Supplementary Table S3). A total of 15,133 KEGG-annotated unigenes were classified into 190 pathways (>10 associated unigenes). Among these pathways, the ten most highly represented were ribosome, carbon metabolism, protein processing in endoplasmic reticulum, oxidative phosphorylation, biosynthesis of amino acids, spliceosome, RNA transport, purine metabolism, peroxisome, and ubiquitin mediated proteolysis (Supplementary Table S4).

Although in most species the male and female genomes differ by a few genes located on sex-specific chromosomes (such as the Y chromosome of mammals), the vast majority of sexually dimorphic traits result from the differential expression of genes that are present in both sexes (Connallon and Knowles, 2005; Rinn and Snyder, 2005; Ellegren and Parsch, 2007), and this is especially true in hymenopteran insects. Because sex determination in hymenopteran species is haplodiploid, females and males are nearly identical genetically (Ellegren and Parsch, 2007). Such DEGs include those that are expressed exclusively in one sex (sex-specific expression) and those that are expressed in both sexes but at a higher level in one sex (sex-biased expression). These sex-biased genes can be further separated into male-biased and female-biased genes, depending on which sex shows higher expression. Genes with equal expression in the two sexes are referred to as unbiased (Ellegren and Parsch, 2007).

Using our assembled transcriptome as a reference, we identified putative genes expressed in males and females using the RPKM method, and genes with at least 2-fold changes and FDR < 0.01 were defined as DEGs. By comparing female and male transcriptomes, 1416 DEGs were found in B. lasus, of which 442 genes were annotated in GO, 420 in COG, 1024 in nr, 613 in Swiss-Prot, and 396 in KEGG (Table 1). Among these DEGs, 986 were up-regulated in females and 430 were up-regulated in males (Supplementary Table S5).

In terms of biological control, parasitoid species have been extensively applied for reducing pest species population sizes (Hassan, 1993; Li, 1994; Terayama, 1999; Zhishan et al., 2003; Parra and Zucchi, 2004; Lim et al., 2006) because parasitoids can propagate on or in other arthropods. The venom of parasitoid wasps, which is injected into a host by females before or at oviposition, is important for the successful development of the progeny. Parasitoid venoms have diverse physiological effects on hosts, including developmental arrest; alteration in growth and physiology; suppression of immune responses; induction of paralysis, oncosis, or apoptosis; and alteration of host behavior (Edwards et al., 2006; Price et al., 2009; Tian et al., 2010; Kryukova et al., 2011). In total, three female-biased genes (c100635.graph_c0, c101314.graph_c0, c101670.graph_c0; Supplementary Table S5) in this study were annotated for venom proteins, which were related to known insect venoms from N. vitripennis and belonged to previously known insect venom families, such as serine proteases (Graaf et al., 2010; Werren et al., 2010). Despite the large diversity of parasitoid wasp species, only a small number of venom proteins have been described from wasps. A wealth of unexplored biomolecules is present in parasitoid venoms; these proteins are of value in basic evolutionary studies, venom biology, host-parasite interactions, and the study of the evolution of life strategies, and they may potentially contain components that could be used in biological control and pharmacology (Moreau and Asgari, 2015).

Transient Receptor Potential channels are cation channels that are mainly considered as unique polymodal cell sensors; TRPs can be subdivided into six main subfamilies: the TRPC (canonical), TRPV (vanilloid), TRPM (melastatin), TRPP (polycystin), TRPML (mucolipin), and TRPA (ankyrin) groups (Gees et al., 2010). Functionally, TRP channels cause cell depolarization when activated, which may trigger many voltage-dependent ion channels. Upon stimulation, Ca²⁺-permeable TRP channels generate changes in the intracellular Ca²⁺ concentration, [Ca²⁺]_i, due to Ca²⁺ entry via the plasma membrane. However, evidence is increasing that TRP channels are also located in intracellular organelles and serve as intracellular Ca²⁺ release channels (Berridge et al., 2000; Bootman et al., 2001; Gees et al., 2010). TRP channels in Drosophila are involved in the perception of sensory signals such as light, temperature, humidity, pheromones, sound, and touch (Lin et al., 2005). In our study, we found 13 TRP channel genes in B. lasus; Nasonia and honey bee contain 12 and 11 genes, respectively, indicating that the number of trp channels seems to be well conserved in Hymenoptera (Werren et al., 2010). Of the TRP channel genes in B. lasus, most belong to two subfamilies: TRPC and TRPA (Table 2).

In Drosophila, TRPC plays an important role in the perception of light signals, i.e., the phototransduction pathway (Leung et al., 2000) (ko: 04745; Supplementary Figure S1-d), which was enriched in B. lasus male adults. In Drosophila, a number of genes in the visual signal transduction pathway have been characterized, with functions including rhodopsin activation, phosphoinoside signaling, and the opening of TRP and TRPL channels (Wolff and Ready, 1993; Zuker, 1996; Leung et al., 2000; Wang and Montell, 2007). Our transcriptional analyses (Figure 2A: FDR < 0.01, log₂ FC = 1.62) and q-PCR results (Figure 2B: t = -3.169, df = 6, p = 0.019), showed that the gene corresponding to trp (c103240.graph_c0) was more highly expressed in B. lasus males, consistent with the phototaxis test. Although both females and males tended to move toward light (Figure 2C: female, Z = -1.34, p < 0.05; male, Z = -1.6, p < 0.05), the tendency to prefer light was significantly influenced by sex in adults (Figure 2C. χ² = 4.17, df = 1, p < 0.05), males more preferring to move to light. This result is similar to the results of research on trp mutants in Drosophila, which had altered phenotypes, including a reduction in light response (Leung et al., 2000; Popescu et al., 2006). Female reduction in light response might be due to their long periods living in the dark to search for hosts and lay offspring into them, as most host species (e.g., pupae of L. dispar or H. cunea) hide in dark environments, such as the litter horizon (Yan et al., 1989; Yang et al., 2001). Surprisingly, five members of the TRPA subfamily, which is involved in sensing environmental temperature, were annotated in our study. Animals must maintain thermal homeostasis and avoid prolonged contact with harmfully hot or cold objects (Caterina, 2007; Karashima et al., 2009). Unlike most parasitoid species, which overwinter in their hosts as eggs or larvae, B. lasus lives through the winter in its adult stage (Yan et al., 1989). Thus, TRPA may be essential for B. lasus adults, allowing them to sense harmful cold during winter. In addition, intraspecific aggregations in B. lasus have been observed in previous research, and an active component that elicited the aggregation response was isolated and identified as 3-hexanone (Mohamed and Coppel, 1987). The effects of aggregation behavior include mating, host attack, defense, and thermoregulation, and in this species, a previous study suggested that aggregation resulted from an increase in reproductive success by increasing the probability of mate location, as well as offering the possibility of mate choice (Mohamed and Coppel, 1987). However, combining the above results, adults may also aggregate at a site for purposes of thermoregulation, especially in winter, in response to cold. Further studies are required to elucidate the nature of this cue.