A total of 416 adults of the ladybird H. axyridis (146 females and 270 males) were collected from an alfalfa field in the campus of Beijing Academy of Agricultural and Forestry Sciences (BAAFS), Beijing, China. The ladybirds were transported to the Institute of Plant and Environment Protection of BAAFS, and maintained in culture cages (45.0 × 50.0 × 50.0 cm) made of aluminum frames and 80 mesh plastic nets. Ladybirds were maintained at a density of 30 pairs per cage. Newly-laid eggs were moved to a new plastic container (30.0 × 40.0 × 30.0 cm) until the first instar larvae emerged. Larvae (30 per cage) were kept in culturing cages with folded paper tape as filler to avoid cannibalism. Newly-eclosed adults were isolated in different plastic containers and separated by sex (the sex of the adults was determined by the color of their mouthpart top: Black/dark = female and white/light = male) and polymorphism (i.e., melanic or succinic phenotype).

Old nymphs of the aphid Myzus persicae were provided daily as food (Ad libitum). Three-day-old virgin succinic females of similar body weight were collected for the experiments. Ladybirds were reared in regulated environments with the following characteristics: T = 25°C, relative humidity = 65%, photoperiod = 16 h light and 8 h dark, illumination intensity = 1,200 lux. All the conditions were managed by an automatic environment regulation system (L-100, Suntech, Beijing, China).

Two libraries were constructed, one for ladybirds that were (and one for ladybirds that were not) exposed to beta-cypermethrin for 24 h. Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad CA, USA) according to the manufacturer's instructions. RNA was then treated with DNase. RNA sample concentration and integrity were determined using a 2100 Bioanalyzer (Agilent Technologies). Poly A-containing mRNAs were enriched using oligo (dT) magnetic beads and then fragmented using RNA Fragmented Reagent. First- and second-strand cDNA was synthesized, purified, and end-repaired. Single nucleotides were added, adapters ligated, and the ligated products purified. Finally, cDNA templates were enriched by PCR amplification. The quality of the library was evaluated and the amount of cDNA quantified using an Agilent 2100 Bioanalyzer and a ABI StepOnePlus Real-time PCR system, respectively. The cDNA was then sequenced at 90 bp using the Illumina HiSeqTM 2000 platform at the Beijing Genomics Institute (BGI, Shenzhen, China).

We first filtered out the sequencing adaptors, the unknown nucleotides larger than 5%, and the low quality reads. The resulting clean reads were assembled using Trinity (Grabherr et al., 2011). The resulting sequences from Trinity assembly were considered as unigenes. For annotation, unigenes were aligned by BLASTx, using an E-value cut-off of 10⁻⁵, against the NCBI non-redundant Swiss-Prot Kyoto Encyclopedia of Genes and Genome (KEGG, http://www.genome.jp/kegg/), and against the Cluster of Orthologous Groups (www.ncbi.nlm.nih.gov/COG) protein database. Gene Ontology of unigenes was analyzed using the Blast2Go software (Conesa et al., 2005). Gene ontology functional classification of all unigenes was performed using the WEGO software (Ye et al., 2006). In addition, unigenes without homology to these databases were forecast for their translation direction and open reading frames (ORF) using the ESTScan software (Iseli et al., 1999).

The relative abundance of transcripts of H. axyridis that were (or were not) exposed to beta-cypermethrin was expressed as fragment per kilobase per million fragment, according to Mortazavi et al. (2008) Differentially-expressed genes between exposed and un-exposed H. axyridis were identified on the basis of the rigorous algorithm false discovery rate ≤ 0.001, and absolute value of log₂ Ratio ≥ 1. Differentially-expressed genes were then subjected to gene ontology enrichment analysis and KEGG pathway enrichment analysis. For gene ontology enrichment analysis, the calculated P-value from the hypergeometric test underwent Bonferroni Correction, and the gene ontology terms with the corrected P-value ≤ 0.05 were significantly enriched in all differentially-expressed genes. For pathway enrichment analysis, pathways with a Q value ≤ 0.05 after multiple testing correction were significantly enriched in all differentially-expressed genes.

The data obtained from the flight, locomotion and respiration tests were analyzed using t-test to separate the means between control and treatment (using ProStat software, Poly Software International).

Figures 2A–D shows the flight performance of succinic female adults of H. axyridis that were (or were not) exposed to a sublethal dose of beta-cypermethrin. The ladybirds that were exposed to beta-cypermethrin showed significantly higher flight frequency (Figure 2A, t = 5.13, df = 90, P < 0.01), longer accumulated flight duration (Figure 2B, t = 3.17, df = 89, P < 0.01), and longer accumulated flight distance (Figure 2C, t = 6.13, df = 90, P < 0.01) as compared with their respective control. By contrast, there were no significant differences of average flight speed between the (Figure 2D, t = 1.94, df = 88, P = 0.57) exposed and un-exposed ladybirds.

Similarly, ladybirds exposed to beta-cypermethrin showed significantly higher walking frequency (Figure 3A, t = 3.1, df = 58, P < 0.01), longer accumulated walking duration (Figure 3B, t = 6.87, df = 78, P < 0.01), longer accumulated walking distance (Figure 3C, t = 8.06, df = 66, P < 0.01), higher average walking speed (Figure 3D, t = 6.67, df = 78, P < 0.01) and a higher frequency of direction change (Figure 3E, t = 2.3, df = 78, P < 0.01) as compared with control.

In addition, ladybirds exposed to beta-cypermethrin had a significantly higher respiratory quotient (Figure 4A, t = 2.23, df = 78, P = 0.014) and respiratory rate (Figure 4B, t = 5.48, df = 78, P < 0.01) as compared with their respective control.

RNA-seq analysis generated approximately a million paired-end reads from the exposed and un-exposed ladybirds. To obtain a comprehensive transcriptional profile of H. axyridis, the total clean reads from gene libraries of exposed and un-exposed ladybirds were combined. De novo assembly produced contigs with a mean length of 59.01 bp. These contigs were further assembled into 65,509 unigenes with an average size of 799.57 bp. There were 7,471 unigenes (14.96%) over 1000 bp in length (Figure S1). The N50 length of the contigs and unigenes was 704 and 795 bp, respectively.

For functional annotation, all unigenes were aligned to the Genbank protein database with a cut-off E-value of 10⁻⁵ using BLASTx. Using this approach, 39.99% of all unigenes (26,020 unigenes) were above the cut-off value. The E-value distribution of the top hits in the non-redundant protein database showed that 41.9% of the unigenes (11,182 unigenes) had significant matches (< 1.0E-45); 58.1% of the matched unigenes had an E-value that ranged from 1.0E-5 to 1.0E-45. As for species distribution, most of the unigene sequences matched best to proteins from the red flour beetle (Tribolium castaneum; 37.44% of the unigenes), followed by the mountain pine beetle (Dendroctonus ponderosae; 5.23% of the unigenes), the pea aphid (Acyrthosiphon pisum; 1.32% of the unigenes) and the Asian citrus psyllid (Diaphorina citri; 1.29% of the unigenes) (Figure S2).

Gene ontology analysis was used to identify the potential function of the predicted proteins. A total of 12,088 unigenes were annotated and assigned to gene ontology terms. The most abundant groups were categorized as follows: Biological process, cellular component, and molecular function. In the biological process category (22,776 unigenes), the most abundant unigene groups were cellular process (5,297 unigenes), and metabolic process (5,100 unigenes). In the cellular component category (10,791 unigenes), the most abundant unigene group was cell part (3,719 unigenes). In the molecular function category (14,781 unigenes), the most abundant unigene groups were catalytic activity (6,588 unigenes) and binding (5,791 unigenes) (Figure S3). These results highlighted the importance of cell communication, metabolic activities, cellular structure, and molecular function in the life-cycle of H. axyridis.

Moreover, to further predict putative protein functions, cluster of orthologous groups analysis were performed. A total 14,258 unigenes were annotated and were classified in 25 categories (Figure S4). Among these 25 cluster of orthologous groups categories, the cluster of “general function prediction only” was the largest group (20.23%; 2,885 unigenes), followed by “replication, recombination, and repair” (8.28%; 1,181 unigenes) and “post-translational modification, protein turnover, chaperones” (7.93%; 1,131 unigenes).

The analysis showed that exposing H. axyridis adults to a sublethal dose of beta-cypermethrin had a significant effect on their transcriptome profile. Most of the differentially-expressed genes were down-regulated (Figure 5). For KEGG enrichment analysis, a total of 415 unigenes were assigned to 103 KEGG pathways, mainly involving metabolism, genetic information processing, environmental information processing, and cellular processes. Among these pathways, 11 pathways were significantly enriched with Q-values < 0.05 (Figure 6). Drug metabolism, starch and sucrose metabolism and neuroactive ligand-receptor interaction were the major enrichment pathways (Figure 7). Among these pathways, two esterase and four glutathione S-transferases were identified. And all of these genes were down-regulation when H. axyridis confronted with beta-cypermethrin after 24 h (Table 1). Taken together, this enrichment analysis indicated that drug metabolic and energy pathways played vital roles in the H. axyridis response to a sublethal dose of beta-cypermethrin.

Flight and locomotion performances are very important for predatory ladybirds. These features are associated with prey hunting. They are also related to migration, colonization, and dispersion in agroecosystems. These features are of great importance when ladybirds face poor environments, such as food shortages, or in high intraspecific density (Zera, 1997). In the present study, we found that exposure to a sublethal dose of beta-cypermethrin induced succinic female adults of H. axyridis to initiate flight more often, and fly for longer times and distances with increased speed. In addition, H. axyridis exposed to beta-cypermethrin showed significantly higher walking frequency and duration, longer walking distance, and higher speed. Pesticides residues may thus stimulate a general motion behavior that includes walking and flying. Sublethal doses of pesticides can impact flight and locomotion of various insects both negatively or positively, depending on the species. The parasitoid wasp Encyrtus saliens increased its walking speed after being exposed to a sublethal dose of malathion (Linn and Roelofs, 1985). The cabbage looper moth Trichoplusia ni presented hyperreflexia after being exposed to a sublethal dose of pesticides and initiated flight more frequently (Hoy and Dahlsten, 1984). By contrast, the green peach aphid M. persicae reduced its willingness to fly after being exposed to a sublethal dose of aldicarb (Boiteau et al., 1985). Sometimes the sublethal effects of pesticides on the activities of insects are indirect. There are numerous studies showing how the flight of the pink bollworm moth Pectinophora gossypiella is influenced by sublethal doses of neurotoxic insecticides (Haynes, 1988). Wiles and Jepson (1994) found that residues of deltamethrin have sublethal effects on predatory ladybirds. Seven spots ladybirds were driven out of wheat agroecosystem containing deltamethrin and migrated to wheat land lacking pesticide. There are two possible explanations for ladybirds to migrate out of agroecosystems containing residual pesticides. On the one hand, the ladybirds may follow the migration of their prey out from the contaminated ecosystem. On the other hand, pesticides may have repellent effects on predators (Desneux et al., 2007). Despite the decrease of the prey densities after the application of pesticides, the acclimation of the ladybird to sublethal doses of pesticides may play also a key role for successful IPM programs.

Any insect activity is related to respiratory and energy metabolism, and is influenced by environmental changes (Wong et al., 2016). The flight of insects is associated with the release of hormones and with the metabolism of lipids, fatty acids, or glycerol (Haunerland, 1997; van der Hourst, 2003). We could not find any reports showing how a sublethal dose of pesticides may alter the energy metabolism benefiting flight and locomotion. Locomotion and flight changes triggered by environmental toxicants could be regarded as external performance modifiers, following responsive physiological paths and complex changes of gene expression. Exploring gene expression changes in response to sublethal doses of pesticides may help us to better understand the physiological responses of the predators to the residual pesticides and to reveal the mechanisms of the toxicant-driven acclimation of ladybirds.

The advent of next-generation sequencing has allowed for rapid sequencing and transcriptomic research. A total of 66,506 unigenes were generated in our study. In addition, the BLASTx analysis of Genbank showed that 39.99% of the unigenes have significant homology to functional genes that encode specific proteins. In consequence, 60.01% of the unigenes had no clear homology to known genes. This low percentage of annotated genes was most likely attributed to the absence of specific information on the genome of H. axyridis. Due to this absence of information, some H. axyridis transcripts, derived from the untranslated regions or from non-conserved domains cannot be annotated.

The red flour beetle, T. castaneum is the only beetle species with a completely sequenced genome (Richards et al., 2008). As expected, T. castaneum is the species that returned most of the BLAST hits (37.44%) against H. axyridis transcripts. Transcriptome sequences were used to carry out research dealing with the effect of insecticides and diet shifts on ladybirds in a previous study (Zhang et al., 2012; Tang et al., 2014; Li et al., 2016). The new transcriptome data of H. axyridis reported here greatly enriches the available genetic information of this ladybird species.

When confronted with insecticides, the detoxification systems of insect are activated to cope with toxic chemicals. Our results demonstrate that a stress through exposure to a sublethal dose of beta-cypermethrin resulted in large variation in the transcriptome profile of ladybirds. Several detoxification pathways were significantly changed. Among these pathways, we detected two esterases and four glutathione S-transferases that were down-regulated. There are several detoxification mechanisms in insects: Cytochrome P450s, glutathione S-transferases, and esterases are recognized as three major detoxifying enzyme systems that are involved in detoxification of various classes of chemicals (Kristensen, 2005; Feyereisen, 2012). Our results suggest that H. axyridis recover from insecticide exposure by eliminating toxic chemicals through detoxifying enzymes. The down-regulation of glutathione S-transferases and esterases could play an important role in this process.

Toxic stress (both from direct spraying of chemical agricultural pesticides and from their residues) may be a strong drive for adaptive evolution in arthropods, which may change their population structures, age, and sexual ratios (Desneux et al., 2007; Guedes et al., 2016). In addition, as shown here, some predators in high trophic levels of the food web could show acclimation by increasing flight or locomotion performances, while also changing their development and reproduction strategies as has been reported elsewhere (Xiao et al., 2016). Transcriptional analysis demonstrated that this kind of plasticity of adaptation may result from the selective expression of functional genes related to energy metabolism and detoxification systems. During adaptive evolution, sublethal impacts of environmental stressors may lead to directional selection, by selecting for the genes that may benefit the population subjected to those specific stressors (Schlichting and Pigliucci, 1998; Li et al., 2014). We must rely on the novel differentially-expressed genes found in this study to explore the possible adaptive evolution of predatory ladybirds facing sublethal doses of toxicants.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer RL and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.