RNA interference-mediated silencing of genes involved in the immune responses of the soybean pod borer Leguminivora glycinivorella (Lepidoptera: Olethreutidae)

Insect rearing

L. glycinivorella eggs collected from a naturally infested soybean field at the experimental station of Northeast Agricultural University in Harbin, China were hatched at 26 °C. The resulting larvae were reared on an artificial diet prepared in our laboratory (Meng et al., 2017a). Adult moths were fed a 5% honey solution, and were allowed to oviposit on young bean pods. The first instar larvae were selected and subjected to artificial diet feeding experiments.

Illumina sequencing

The T3 dsSpbP0 (double-stranded SBP ribosomal protein P0 RNAi transgenic soybean line) and wild-type ‘DN50’ plants, provided by the Key Laboratory of Soybean Biology of the Chinese Education Ministry, Harbin, China, were grown in a greenhouse at 24 ± 1 °C with 60% relative humidity under a 16-h light/8-h dark cycle (Meng et al., 2017b). At the R5 soybean stage (fully developed pods), three replicates of 50 first-instar larvae were reared on soybean pods of DN50 and T3 dsSpbP0 plants. All larvae were collected after three d and used to construct a cDNA library. cDNA library preparation and Illumina sequencing were conducted by the Biomarker Technology Company (Beijing, China). Briefly, total RNA was extracted from six pooled larvae using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The first-strand cDNA was synthesised using random hexamer-primers from purified Poly (A) mRNA. Second-strand cDNA was synthesised using DNA polymerase I and RNaseH, and then purified using a QiaQuick PCR extraction kit (Qiagen, Hilden, Germany). cDNA fragments of a suitable length (300–500 bp) were obtained by agarose gel electrophoresis and amplified by PCR to construct the final cDNA libraries using NEB Next Ultra RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA). The cDNA library was sequenced on the Illumina HiSeq 2000 system (Illumina, San Diego, CA, USA). The unigenes from six samples were combined to create the SPB unigene database (Chen et al., 2014; Meng et al., 2017c). All raw transcriptome data have been deposited in the NIH Short Read Archive (accession numbers SRR5985984, SRR5985985, SRR5985986, SRR5985987, SRR5985988, SRR5985989).

Identification of immune-related genes

A list of immune-related genes was compiled based on the available relevant literature (Table S1), and homologous B. mori, D. melanogaster, Danaus plexippus and Papilio polytes genes in the GenBank database were identified (Guan & Mariuzza, 2007; Xu et al., 2012). The tBLASTn algorithm-based tool was used to complete sequence similarity searches of the SPB transcriptome database (Boratyn et al., 2013). The 41 SPB immune-related genes were submitted to NCBI GenBank and their Accession numbers are shown in Table S1.

Phylogenetic and domain analyses

Amino acid sequences were aligned with the Multiple Alignment program clustal omega (https://www.ebi.ac.uk/Tools/msa/clustalo/), and the phylogenetic tree was constructed in MEGA 5 based on the neighbour-joining method with 1,000 bootstrap replicates (Tamura et al., 2011; Sievers et al., 2011). The architecture of the protein domains was analysed using the SMART program (http://smart.embl-heidelberg.de/).

dsRNA synthesis

We synthesised dsRNAs using the T7 RiboMAX Express Large Scale RNA Production System (Promega, Madison, WI, USA). The T7 RNA polymerase promoter sequence was added to each end of the DNA templates during PCR amplifications. The primers containing the T7 RNA polymerase promoter were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (Table S2). For the negative control, the green fluorescent protein (GFP) gene was amplified from the PCAMBIA1302 vector as a template for GFP dsRNA synthesis. The template DNA and single-stranded RNA were eliminated from the transcription reaction by DNase I and RNase A, respectively. The prepared dsRNAs were purified by a phenol/chloroform extraction followed by an ammonium acetate precipitation. The dsRNAs were ultimately suspended in ultrapure water and quantified using a Nano Drop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA).

SPB feeding bioassay

The first instar SPB larvae were fed an artificial diet containing dsRNA (10 µg/g) for specific target genes as described by Meng et al. (2017a). Control larvae were treated with the same concentration of GFP dsRNA. The feeding bioassay was completed in triplicate with 50 larvae per treatment or control. Three biological replicates were used for each treatment. The larvae were reared for 15 d at 26 °C under a 16-h light/8-h dark cycle with 65% relative humidity. The dsRNA-supplemented artificial diet was refreshed every three d. Body weight, mortality, and phenotypic abnormalities were recorded every three d.

Quantitative real-time PCR (qRT-PCR)

For every treatment, two surviving larvae were randomly collected at each time point from all of the biological replicates from zero to 15 d after larvae were fed the artificial diet containing the target gene’s dsRNA, frozen in liquid nitrogen and stored at −80 °C. Total RNA was extracted from two pooled larvae using the RNApure Tissue Kit (DNase I) (CWBIO, Beijing, China). Additionally, primer sets were synthesised (Table S2). The extracted RNA samples were treated with DNase I (Invitrogen, Carlsbad, CA, USA) to remove any contaminating genomic DNA. They were then used as the template for first-strand cDNA synthesis with the TIANScript RT Kit (Tiangen, Beijing, China). The qRT-PCR analysis was completed using the SYBR Green kit (Bio-Rad, Hercules, CA, USA) and a Roche LightCycler® 480 real-time PCR system (Roche, Basel, Switzerland). Each reaction used cDNA corresponding to 50 ng of total RNA and each primer at a final concentration of 100 nm in 20 µl reaction. Controls included non-RT controls (50 ng total RNA without reverse transcription was used to detect genomic DNA contamination) and non-template controls (water template).The qRT-PCR conditions were as follows: 95 °C for 5 min; 40 cycles at 95 °C for 30 s, 60 °C for 15 s, and 72 °C for 45 s; 95 °C for 1 min and 55 °C for 1 min. At the end of each qRT-PCR experiment, a melt curve was generated to check for primer-dimer formation. The efficiencies of the qRT-PCR primer pairs were greater than 90% (Table S2). The qRT-PCR analysis contained three biological replicates, each having three technical replicates (Meng et al., 2017c). Each target gene’s qRT-PCR products were sequenced to confirm their identities.

Relative expression levels were calculated by the following formula using LightCycler 480 Software v1.5.0 (Roche, Basel, Switzerland): R = 2^{−(ΔCt sample−ΔCt calibrator)}, where R represents the relative expression level, ΔCt sample is the average difference between the Ct of the gene and that of SBPß-actin in the experimental sample (Meng et al., 2017b), and ΔCt calibrator is the average difference between the Ct of the gene and the e SBPß-actin in the calibrator. A representative sample was set as the calibrator (Bustin et al., 2009). The differences in Ct values between technical replicates was less than 0.5. The relative expression level of 0d for each target gene was set as a benchmark. The relative expression level of other time point for each target gene was used to compare with 0d value.

Statistical analyses

All of the data in this study are presented as mean ± SE. Significant differences were determined by one-way analysis of variance followed by least significant difference tests for mean comparisons. The statistical analysis was performed with SAS 9.21 software (SAS Institute, Cary, NC, USA). P-values were corrected by Bonferroni.

Identification of iLeguminivora glycinivorella immune-related genes

Amino acid sequences encoded by D. melanogaster, B. mori, and Manduca sexta (Lepidoptera) immune-related genes were used to search SPB transcriptome sequences. The 41 identified putative SPB immune-related genes were functionally classified into three groups, microbial recognition, immune signalling, and immune effector molecules (Table S1).

Microbial recognition molecules

The PGRPs recognize conserved molecular patterns present in pathogens, but absent in the host, including PGNs, which are essential cell wall components of almost all bacteria. The PGRPs are encoded by a highly conserved gene family in insects, and are generally classified into two types (short and long) (Dziarski & Gupta, 2006; Yang et al., 2017). We identified eight SPB PGRPs, four short and four long types, with similarities to D. melanogaster PGRP-SC, PGRP-SD, and PGRP-LB (Fig. 1). Five (LgPGRP-SC1a, LgPGRP-SC1b, LgPGRP-SD1b, LgPGRP-LB1, and LgPGRP-LB2b) of the eight identified PGRPs were predicted to be secreted proteins, based on the presence of putative signal peptides, that function as amidases. LgPGRP-LB lacks a putative signal peptide, but contains a transmembrane region and amidase domain, suggesting that it serves as a transmembrane PGN receptor. In contrast, LgPGRP-SD1a and LgPGRP-LB2a carry only the PGRP domain, implying that they are intracellular proteins (Table S1).

The gram-negative bacteria-binding proteins (GNBPs) and β-1,3-glucan recognition proteins (βGRPs) belong to a subfamily of pattern recognition receptors and have strong affinities for the lipopolysaccharide of gram-negative bacteria and the β-1,3-glucan of fungi, but not for the PGN of gram-positive bacteria. Of the three GNBPs produced by D. melanogaster (GNBP1, GNBP2, and GNBP3), GNBP1 interacts with PGRP-SA to form a hydrolytic complex that activates the Toll pathway in response to gram-positive bacteria, while GNBP3 is required for detecting fungi and activating the Toll pathway (Hughes, 2012; Rao et al., 2017). We identified one GNBP gene in the SPB transcriptome datasets. The neighbour-joining phylogenetic analysis indicated that LgGNBP3 is an ortholog of DmGNBP3 (Fig. S1). A comparison between the deduced amino acid sequences and the D. melanogaster GNBP sequences indicated that LgGNBP3 contains a putative N-terminal β-1,3-glucan-recognition domain (CBM39) and a C-terminal glucanase-like domain (glycosyl hydrolase family 16), suggesting that LgGNBP3 may bind to fungal β-1,3-glucan.

Immune signalling molecules

After specific ligands are detected, microbial recognition molecules activate or modulate various immune response pathways. Here, we identified genes associated with the Toll and IMD pathways, which are the major signalling pathways that mediate the innate immunity of insects. The Toll pathway regulates the production of antimicrobial peptides in response to infections by fungi or gram-positive bacteria with lysine-type PGNs in their cell walls (Roh et al., 2009). The Toll receptor, which is responsible for the signal transduction associated with the Toll pathway, is vital for insect innate immune responses and embryo development (Takeda & Akira, 2004; Benton et al., 2016). In this study, we identified 10 genes encoding Toll receptors in the SPB transcriptome datasets. The TIR domain is highly conserved in insect Toll families. To investigate the orthologous relationships among these genes, we constructed a phylogenetic tree based on an alignment of the TIR domains from all SPB and D. melanogaster Toll proteins. The Toll receptors analysed in this study formed six major clusters, namely Toll-3, Toll-4, Toll-5, Toll-6, Toll-7 and Toll-8 (Fig. 2A). Based on the phylogenetic tree, the L. glycinivorella Toll genes were designated as LgTLR-3, LgTLR-4, LgToll-5-1a, LgToll-5-1b, LgToll-6, TLR-6-1, TLR-6-2, LgToll-7, LgTLR-7 and LgToll-8. All 10 predicted proteins contain an extracellular LRR domain as well as transmembrane and cytoplasmic TIR domains (Fig. 2B). We also identified sequences matching the intracellular components, ECSIT and Tollip, which affect the Toll signalling pathway (Table S1).

The IMD pathway is mainly activated by gram-negative bacterial infections. Additionally, IMD signal transduction is reportedly mediated by IMD, fas-associated death domain protein (FADD), death-related ced-3/Nedd2-like caspase (Dredd), inhibitor of apoptosis protein 2 (IAP2), transforming growth factor β-activated kinase (TAK1), TAK1-binding 2 (Tab2), ubiquitin conjugating 13 (Ubc13), and an inhibitor of nuclease factor B kinase subunits b and g (IKKb and IKKg) (Bao et al., 2013; Myllymäki, Valanne & Rämet, 2014). Of these,we only identified sequences that were homologous to FADD, Dredd, and IAP2 (Table S1).

Immune-related effector genes

The PGRPs and βGRPs detect PGNs and β-1,3-glucans, which activate a clip-domain serine protease (CLIP) cascade that converts prophenoloxidase to active phenoloxidase, leading to the melanisation responses involved in eliminating pathogens (Monwan, Amparyup & Tassanakajon, 2017; Li et al., 2016). We identified two CLIP genes (LgSnake-1 and LgSnake-2) in the SPB transcriptome datasets. The deduced amino acid sequences each contain a clip domain at the N-terminus and a serine protease domain at the C-terminus (Fig. S2). Serine protease inhibitors (i.e., serpins) negatively regulate prophenoloxidase activation, which prevents the excessive activation of the CLIP cascade. In this study, we identified three serpin genes, Lgserpin1, Lgserpin2, and Lgserpin3, in the SPB transcriptome datasets. Their deduced amino acid sequences each contained a putative signal peptide sequence and a core serpin domain, suggesting that they are secreted proteins (Table S1).

Immune response effector genes

Antibacterial peptides are immune response effectors that are induced by immune challenges and are important for defence responses against insects. Diverse antibacterial peptide genes have been identified in many insect species, including genes encoding defensins, reeler, and lysozyme (Imler & Bulet, 2005; Bao et al., 2011). In this study, we identified two defensin genes, Lgdefensin1 and Lgdefensin2, in the SPB transcriptome datasets. The encoded amino acid sequences each consisted of a putative signal peptide sequence and a core Knot1 domain (Table S1). We also identified seven chicken-type (C-type) lysozymes and two invertebrate-type (I-type) lysozymes in the L SPB transcriptome. The C-type lysozymes are bacteriolytic enzymes that hydrolyse the β (1–4) bonds between N-acetylglucosamine and N-acetylmuramic acid in the PGN of prokaryotic cell walls. The predicted SPB C-type proteins, with the exception of the C-type 3 protein, each include an N-terminal signal peptide sequence (Table S1). Additionally, we detected eight conserved cysteine residues in the L. glycinivorella C-type lysozymes (Fig. S3A) as well as 12 conserved cysteine residues in the deduced SPB I-type lysozyme sequences (Fig. S3B). These cysteine residues may form intramolecular disulphide bonds to enhance stability and resistance against proteolytic degradation.

Potential RNAi targets identified in an artificial feeding assay and effects of double-stranded RNA on soybean pod borer development and mortality

In total, 11 genes representing the immune-related SPB genes were selected and analysed to identify potential new RNAi targets useful for controlling the SPB. We synthesised the corresponding dsRNAs in vitro and mixed them in an artificial diet. The mortality rates three d after larvae were fed artificial diets, which contained containing 10 µg/g dsRNA for LgPGRP-LB, LgPGRP-LB2b, LgToll-5-1a, LgToll-5-1b, LgTLR-7, LgSerpin2 or LgChaoptin were 37%–92%. These mortality rates were significantly greater than those of control larvae treated with GFP dsRNA. Additionally, the final mortality rates of the larvae fed LgPGRP-LB, LgPGRP-LB2b, LgToll-5-1a, LgSerpin2 and LgChaoptin dsRNA were even higher at 15 d. In contrast, the artificial diets containing dsRNA independently targeting LgPGRP-LB2a, LgTLR-3 and Lgitype-1 did not have any statistically significant effects on larval mortality (Fig. 3). Moreover, two main phenotypic differences were observed among the surviving larvae after 15 d of feeding. First, the weights of the larvae fed LgToll-5-1b and LgItype-1 dsRNA increased more gradually than those of larvae fed GFP dsRNA, and was ultimately lower after 15 d of feeding (Fig. 4). Additionally, the LgToll-5-1b RNAi treatment resulted in partly black cuticles. Second, larvae fed dsRNA targeting LgPGRP-LB or LgPGRP-LB2a underwent early pupation, with pupation rates of 25 and 50%, respectively. The remaining larvae developed abnormally with stunted and twisted bodies (Fig. 5).

To investigate how larval mortality and abnormal development are correlated with the relative expression levels of specific target genes, we performed qRT-PCR using total RNA extracted from the surviving larvae at different time points after feeding on artificial diets. The expression levels of the genes, except for those of LgPGRP-LB2a, LgTLR-6-1 and LgChaoptin decreased in larvae three d after being treated with the respective dsRNAs, while the expression levels of all of the genes decreased significantly after six or nine d (Fig. 6). Thus, the increased mortality rates and abnormal development of larvae fed dsRNA may result from the down-regulated expression of specific target genes. Moreover, unigenes LgToll-5a and LgPGRP-LB2a may represent good RNAi targets for controlling the SPB.