Introduction
MicroRNAs (miRNAs) are single-stranded, endogenous, small non-coding RNAs, with approximately 22 nucleotides in length (nt) (Lai et al., Reference Lai, Wiel and Rubin2004; Vaucheret et al., Reference Vaucheret, Vazquez, Crété and Bartel2004; Cai et al., Reference Cai, Yu, Hu and Yu2009). These small molecules derive from the primary miRNA (pri-miRNA) transcribed by RNA polymerase II Drosha (Lee et al., Reference Lee, Ahn, Han, Choi, Kim, Yim, Lee, Provost, Rådmark, Kim and Kim2003; Denli et al., Reference Denli, Tops, Plasterk, Ketting and Hannon2004). The pri-miRNAs undergo a series of processing and sorting events in the nucleus to form precursor miRNA (pre-miRNA) (Lucas et al., Reference Lucas, Zhao, Liu and Raikhel2015). Once the pre-miRNAs are transported to the cytoplasm, they can be cleaved by RNase III endonuclease, Dicer1, to produce a 20–25 nt length miRNA–miRNA* duplex (Llave et al., Reference Llave, Xie, Kasschau and Carrington2002). The duplex is loaded into Argonaute 1 protein (Ago1) which can form an RNA-induced silencing complex (RISC), and one strand is retained for functioning (Filipowicz et al., Reference Filipowicz, Bhattacharyya and Sonenberg2008). The RISC loaded with miRNA could regulate transcript levels of protein-coding genes by pairing to the specific site of messenger RNA (mRNA) (Fang and Rajewsky, Reference Fang and Rajewsky2011; Lucas et al., Reference Lucas, Zhao, Liu and Raikhel2015). In insects, numerous functions of miRNAs have been proven to be involved in the development and reproduction, immunity to entomopathogen and susceptibility to insecticides (Yang et al., Reference Yang, Wei, Jiang, Wang, Guo, He and Kang2014; Wei et al., Reference Wei, Zheng, Peng, Pan, Xi, Chen, Zhang, Yang, Gao and Shang2016; Ma et al., Reference Ma, Liu, Zhao, Yang, Chen, Li and Lu2020). For instance, Helicoverpa armigera miR-2055 directly regulates lipid metabolism via fatty acid synthase, then indirectly affects development and reproduction (Cheng et al., Reference Cheng, Lu, Guo, Lin, Jin, Zhang and Zou2022). In addition, Spodoptera exigua miR-998-3p causes the up-expression of ATP-binding cassette transporter proteins ABCC2 to decrease sensibility to Cry1Ac (Zhu et al., Reference Zhu, Sun, Nie, Liang and Gao2020).
As one kind of entomopathogen, Bacillus thuringiensis, a Gram-positive bacterium, could produce insecticidal proteins (i.e. Cry, Cyt and Vip proteins) and has been widely used to manage pests (Crickmore et al., Reference Crickmore, Berry, Panneerselvam, Mishra, Connor and Bonning2021; Tabashnik et al., Reference Tabashnik, Liesner, Ellsworth, Unnithan, Fabrick, Naranjo, Li, Dennehy, Antilla, Staten and Carrière2021; Ji et al., Reference Ji, Gao, Zhao, Wang, Zhang, Wu, Xie, Shi and Guo2024). Currently, studies about B. thuringiensis mainly focused on insecticidal toxicology. There are two proposed models to elucidate Bt insecticidal mechanisms: the first is that Bt causes an osmotic imbalance in response to the formation of pores in a midgut epithelial cell membrane, and the second is that it causes an opening of ion channels that activate the process of cell death, which in turn leads to insect death (Sanahuja et al., Reference Sanahuja, Banakar, Twyman, Capell and Christou2011; Melo et al., Reference Melo, Soccol and Soccol2016; Bel et al., Reference Bel, Ferré and Hernández-Martínez2020). In addition to the direct effect, B. thuringiensis can also affect the growth, development and immune response of insects (Hussein et al., Reference Hussein, Habuštová and Sehnal2005; Zhang et al., Reference Zhang, Ma, Wan, Mu and Li2013b; Grizanova et al., Reference Grizanova, Dubovskiy, Whitten and Glupov2014). Previous study showed that the weight of Mythimna separata decreased after feeding corn leaves covered with Cry1Ac or Cry2Ab protein (Wang et al., Reference Wang, Lu, Li, Li, Dong and Liu2018). The adult emergence of Acanthoscelides obtectus that fed on the diet with Cry1Ia, Cry 7Ab and Cry23/37 proteins was lower, compared with the control groups (Rodríguez-González et al., Reference Rodríguez-González, Porteous-Álvarez, Del Val, Casquero and Escriche2020). Recent evidence suggests that miRNAs can be involved in the immune response of insects caused by B. thuringiensis infection. The question for whether miRNAs can be involved in the other detrimental effects caused by B. thuringiensis, and what kinds of miRNAs can regulate these effects still remains.
The beet armyworm, S. exigua (Hübner), is an important pest worldwide (Greenberg et al., Reference Greenberg, Sappington, Legaspi, Liu and Sétamou2001; Feng et al., Reference Feng, Wu, Cheng and Guo2003; Rabelo et al., Reference Rabelo, Santos and Paula-Moraes2022). It could feed on vegetables, flower crops and other agricultural crops (Moulton et al., Reference Moulton, Pepper and Dennehy2000; Maharjan et al., Reference Maharjan, Ahn and Yi2022). Traditional chemical pesticides led to insect resistance and environmental pollution (Hafeez et al., Reference Hafeez, Ullah, Khan, Li, Zhang, Shah, Imran, Assiri, Fernández-Grandon, Desneux, Rehman, Fahad and Lu2022; Rabelo et al., Reference Rabelo, Santos and Paula-Moraes2022). The application of B. thuringiensis can effectively manage S. exigua. However, our understanding of S. exigua miRNAs in mediating B. thuringiensis infection is not only limited, but also the characteristics, expression and functions of miRNA core genes of S. exigua infected by B. thuringiensis are still unknown.
In this study, we constructed small RNA (sRNA) libraries of S. exigua infested by B. thuringiensis GS57 strain to confirm whether and what kinds of miRNAs are involved in regulating a variety of insect biological processes caused by B. thuringiensis infection using Illumina sequencing. To explore miRNA functions, the differentially expressed miRNAs (DE miRNAs) and potential target genes of those were analysed. The identification and characteristics of miRNA core genes, SeDicer1 and SeAgo1 were constructed. We further confirmed that the B. thuringiensis infection affected the expression of SeDicer1 and SeAgo1. The knockdown of SeDicer1 directly decreased the expression level of miRNAs, suggesting the essentiality of SeDicer1 in miRNA regulation, and the survival rate of S. exigua infected by B. thuringiensis GS57 decreased after injecting dsSeDicer1. These data provide insight into understanding the relationship between S. exigua and B. thuringiensis via miRNA, which can help pest management strategies in the future.
Materials and methods
Insect larvae and B. thuringiensis strain
A S. exigua colony was provided by the Jilin Academy of Agricultural Sciences, China. Larvae were reared with an artificial diet which was followed by the method described previously, under controlled conditions of 27 ± 1°C and 16L: 8D h (Ren et al., Reference Ren, Chen, Zhang, Ma, Cui, Han, Mu and Li2013). The colony was maintained for more than 20 generations without exposure to any insecticides and B. thuringiensis strains in the laboratory.
The B. thuringiensis GS57 strain was isolated from soil and maintained in our laboratory. This strain was isolated by using the temperature screening method (Su et al., Reference Su, Zhang, Yuan, Li, Ma and Tan2007). Briefly, 0.1 g of soil sample was put into a test tube with 10 ml sterilised water and glass beads. Then, the tube was shaken with 200 rpm for 20 min, and water bathed at 75°C for 20 min to ensure the inactivation of non-bacillus bacteria. After standing for 1 min, 100 μl bacterial solution which was diluted to the concentration of 10−2, 10−3 and 10−4 was coated on the 1/2 LB solid medium, respectively. Then, the solid mediums were cultured at 30°C for 3 days. Suspected B. thuringiensis colonies were selected and examined with carbolic acid red staining. The insect toxicity of B. thuringiensis GS57 against S. exigua larvae has been shown in our previous studies (Li et al., Reference Li, Zhao, Wu, Ji, Liu, Guo, Guo and Bi2022b). Bacillus thuringiensis GS57 was inoculated in 1/2 LB medium and incubated for 46 h at 30°C until 70–90% of crystals were released (Li et al., Reference Li, Zhao, Wu, Ji, Liu, Guo, Guo and Bi2022b). The bacteria and crystals were centrifuged and subsequently re-suspended in sterile double-distilled water (ddH2O) with the final concentration of 10 mg ml−1 for assay.
RNA sample, small RNA library construction and sequencing
The fourth instar larvae of S. exigua were reared on the artificial diet covered with B. thuringiensis GS57 of 10 mg ml−1. Spodoptera exigua larvae treated with B. thuringiensis GS57 (Bt) at 24, 48, 72 and 96 h and sterile ddH2O (CK) were randomly sampled for extracting RNA. All experiments were performed by three biological replicates of five larvae.
According to guidelines of the manufacturer's instructions, total RNA was isolated from Bt and CK using Trizol reagent (Invitrogen, Carlsbad, CA, USA). Firstly, samples were immersed and homogenised in liquid nitrogen. Then, 50–100 mg were used for RNA isolation. A total of 1 ml Trizol reagent was added to the tube, and samples were homogenised using a power homogeniser. The homogenised sample was incubated for 5 min at 25°C. After adding 0.2 ml chloroform, the tube was vigorously shaken for 15 s and incubated for 3 min at 25°C. Followed by centrifuging at 12,000 × g for 15 min at 4°C, the aqueous phase of the sample was transferred to a new tube, and then 0.5 ml isopropanol was added into the aqueous phase and homogenised, incubated at 25°C for 10 min and then centrifuged at 12,000 × g for 10 min at 4°C. The supernatant was removed from the tube, leaving the RNA pellet, washed the pellet with 1 ml of 75% ethanol and then centrifuged the tube at 7500 × g for 5 min at 4°C. The wash solution was then removed. There were two repetitions of the pellet washing. To completely remove the ethanol, the RNA pellet was dried in the air for 10 min. Lastly, the RNA was resuspended in RNase-free water for downstream application. The concentration and purity of RNA were determined using the NanoDrop 2000 (Thermo Fisher Scientific, USA) and Agient2100, LabChip GX (PerkinElmer, USA), respectively.
The sRNA libraries of Bt and CK (treated with 72 h) were constructed, amplified, sequenced and analysed by Illumina at BioMarker Biotechnology Co., LTD (BioMarker, China) using VAHTS Small RNA Library Prep Kit for Illumina (Vazyme, China). Firstly, the 3’ SR and 5’ SR adaptors were ligated to RNA, and then transcription was reversed to synthesise first chain. Lastly, PCR amplification and size selection were conducted. PAGE gel was used for electrophoresis fragment screening purposes, and rubber cutting recycling was used as the pieces get sRNA libraries. At last, PCR products were purified (AMPure XP system) and library quality was assessed.
For sequencing, the clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v4-cBot-HS (Illumina, USA). Based on cluster generation, the library preparations were sequenced on an Illumina platform and single-end reads were generated.
Bioinformatics analysis of small RNA sequences
Raw reads were firstly processed through in-house Perl scripts. In this step, clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data; and reads were trimmed and cleaned by removing the sequences smaller than 18 nt or longer than 30 nt which is consistent with length distribution of sRNA. At the same time, Quality Score 20, Quality Score 30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.
Using Bowtie tools soft (v1.0.0) (parameter: v, 0), the clean reads were aligned with Silva (http://www.arb-silva.de/), GtRNAdb (http://lowelab.ucsc.edu/GtRNAdb/), Rfam (http://rfam.xfam.org/) and Repbase (http://www.girinst.org/repbase/) database to filter rRNA, tRNA, snRNA, snoRNA and other ncRNA and repeats. The remaining reads were used to predict known and novel miRNAs by comparing with the S. exigua genome sequence (GCA_022117675.1) (Simon et al., Reference Simon, Breeschoten, Jansen, Dirks, Schranz and Ros2021). For known miRNA prediction, reads matched to the S. exigua genome were compared with the sequences of known miRNAs in miRBase (v22) with one mismatch allowed (Ambros et al., Reference Ambros, Bartel, Bartel, Burge, Carrington, Chen, Dreyfuss, Eddy, Griffiths-Jones, Marshall, Matzke, Ruvkun and Tuschl2003). Using miRDeep2 soft (v2.0.5) (parameter: g, −1; b, 0), the potential precursor sequences were obtained. The distribution information of reads on the precursor sequences, characteristics of miRNA production (mature sequence, star sequence, loop structure) and energy information of precursor structure were used to predict the novel miRNAs. Randfold tools soft (parameter: s, 99; default) was used for the energy of pre-miRNA structure and novel miRNA secondary structure prediction.
The expression level of miRNAs was estimated by transcripts per million (TPM) (Love et al., Reference Love, Huber and Anders2014). The analysis of DE miRNAs between Bt and CK were performed using the DESeq2 R package (v1.10.1). The package DESeq2 provides statistical routines for determining differential expression in digital miRNA expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. miRNA with |log2(FC)|≥0.58; P≤0.05 was assigned as differential expression.
Moreover, miRanda (v3.3a) (parameter: sc, 50; en, −20; scale, 4; go, −2.0; ge, −8.0) and targetscan (v5.0) (parameter: default) were used to predict and analyse potential target genes of DE miRNAs (Krüger and Rehmsmeier, Reference Krüger and Rehmsmeier2006; Kuhn et al., Reference Kuhn, Martin, Feldman, Terry, Nuovo and Elton2008). Gene Ontology (GO) database (http://www.geneontology.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses (http://www.genome.jp/kegg/) were conducted to predict target genes functions.
Real-time quantitative PCR (RT-qPCR) analysis of mRNA and miRNA
The first-strand cDNA synthesis of miRNA and mRNA was conducted using a miRcute Plus miRNA First Strand Synthesis CDNA Kit (TIANGEN, China) and PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, Japan), according to the manufacturer's instructions. RT-qPCR of miRNA and mRNA was carried out using a miRcute Plus miRNA qPCR Kit (SYBR Green) (TIANGEN), and SYBR Premix Ex Taq (TaKaRa), following the instructions of the manufacturer. RT-qPCR amplification of mRNA was conducted in a 20 μl reaction volume consisting 10 μl of TB Green Premix Ex Taq 2, 1 μl each of forward and reverse primer, 1 μl of 10× diluted cDNA template and 7 μl ddH2O. A three-step PCR was employed for amplification, with cycling parameters as follows: 1 cycle of 95°C for 30 s, followed by 40 cycles of 95°C for 10 s, 58°C for 30 s and 72°C for 30 s. RT-qPCR amplification of miRNA was conducted in a 20 μl reaction volume consisting 10 μl of miRcute Plus miRNA PreMix, 0.4 μl each of forward and universal reverse primer, 1 μl of 50 × diluted cDNA template and 8.2 μl ddH2O. A two-step PCR was employed for amplification, with cycling parameters as follows: 1 cycle of 95°C for 15 min, followed by 40 cycles of 95°C for 10 s, 58°C for 30 s and 72°C for 30 s. The RT-qPCR was performed on a CFX96 System (Bio-Rad, USA) with three biological replicates. The previous studies have proven that the expression of β-actin is stable under B. thuringiensis infection (Park and Kim, Reference Park and Kim2013), and small nuclear RNA U6 is stable under different conditions (Zhu et al., Reference Zhu, Sun, Nie, Liang and Gao2020; Liu et al., Reference Liu, Zhang, Fang, Wu and Weng2022). Thus, small nuclear RNA U6 and β-actin genes were used as reference gene for normalising the expression level of miRNA and mRNA, respectively. Primers for RT-qPCR were shown in table S1.
Characterisation of miRNA core genes
The previous studies suggested that Dicer1 and Ago1 were core genes involved in the processing, synthesis and function of miRNA (Rahimpour et al., Reference Rahimpour, Moharramipour, Asgari and Mehrabadi2019; Jouravleva et al., Reference Jouravleva, Golovenko, Demo, Dutcher, Hall, Zamore and Korostelev2022; Lee et al., Reference Lee, Lee, Kim, Kim and Roh2023). To determine the effect of B. thuringiensis infection on the expression levels of miRNA core genes, the genome of S. exigua (GCA_022117675.1) was obtained from the National Centre for Biotechnology Information (NCBI). To obtain protein sequence of S. exigua Dicer1 and Ago1, the sequences of Drosophila melanogaster (GenBank: NP_524453.1), Spodoptera litura (GenBank: XP_022832341.1), Manduca sexta (GenBank: XP_037296641.1), Tribolium castaneum (GenBank: XP_008199045.1), Amyelois transitella (GenBank: XP_013188945.1) and Papilio xuthus (GenBank: KPJ05873.1) were used as query Dicer1 sequences, and Ago-1, Bombyx mori (GenBank: NP_001095931.1), D. melanogaster (GenBank: NP_001246314.1), Samia ricini (GenBank: AID68365.1), Mayetiola destructor (GenBank: AFX89034.1), Nilaparvata lugens (GenBank: AGH30326.1), S. litura (GenBank: AHC98009.1), Blattella germanica (GenBank: CCV01212.1) and T. castaneum (GenBank: KYB26000.1) were used as query Dicer1 sequences. Basic Local Alignment Search Tool (BLAST) (e-value <10−5) (http://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to search homologous sequences of miRNA core genes of S. exigua, SeDicer1 and SeAgo1 (Rahimpour et al., Reference Rahimpour, Moharramipour, Asgari and Mehrabadi2019). The alignments of SeDicer1 and SeAgo1 were performed using MAFFT (v7.0) according to the E-INS-i iterative refinement methods (https://mafft.cbrc.jp/alignment/server/). Protein domains were analysed by SMART (http://smart.embl-heidelberg.de/) and Pfam (http://pfam.xfam.org/). According to the Poisson model, the phylogenetic trees were conducted by the MEGA X using the neighbour-joining (NJ) method with 1000 times bootstrap sampling.
RNA interference
In this study, double-stranded RNA (dsRNA) was used to silence the SeDicer1 gene. The dsGFP (GenBank: KJ668651.1) was synthesised to avoid the effect on RNAi by injection of dsRNA. Primers for dsRNA were showed in table S1. The dsDicer1 and dsGFP were synthesised according to the instruction of T7 RiboMAXTM Express RNAi System (Promega, USA). The early fourth instar larvae were starved for 2 h and frozen on ice for 10 min before injection (Ji et al., Reference Ji, Gao, Zhao, Wang, Zhang, Wu, Xie, Shi and Guo2024). Using a manual 1701RN-microinjector (Hamilton, Romania), 2 μl (10 μg) of dsRNA was injected at the third abdomen leg into the fourth instar S. exigua larvae cavity followed by rearing on artificial diets covered with B. thuringiensis GS57 of 10 mg ml−1. After 24 h infection, 2 μl (10 μg) of dsRNA was second injected into the S. exigua larvae to keep the silence of SeDicer1. During injection, no fluid outflow was considered as the basic criterion. Finally, S. exigua larvae treated with dsRNA were collected at 72 h post infection (hpi) with B. thuringiensis GS57 of 10 mg ml−1 (after the first injection). The decrease of expression of gene was recognised as the gene silencing (Zhang et al., Reference Zhang, Lu and Zhou2013a; Ji et al., Reference Ji, Gao, Zhao, Wang, Zhang, Wu, Xie, Shi and Guo2024). Eighty larvae injected with dsRNA were reared on the artificial diet covered with B. thuringiensis GS57 of 10 mg ml−1 to determine the effect of dsRNA and B. thuringiensis GS57 on survival rate.
Data analysis
CT values were the average of the three technical replicates and three biological replicates. The data of relative expression level of miRNA and mRNA have been calculated using 2–ΔΔCt (Livak and Schmittgen, Reference Livak and Schmittgen2001; Pfaffl, Reference Pfaffl2001). Between Bt and CK groups, the difference of expression level of SeDicer1, SeAgo1, potential target genes mRNA and DE miRNA were analysed using Student's t test. The difference of expression level of SeDicer1 and DE miRNAs between dsGFP and dsSeDicer1 groups was also analysed using Student's t test. These data were shown as the mean ± SE. These statistical analyses were conducted using the SAS (v8.1) (SAS Institute, Cary, NC, USA) and plotted with GraphPad Prism 8 software. The difference in the survival rate of S. exigua between dsGFP and dsSeDicer1 was analysed using Logrank Mantel-cox test. This statistical analysis was conducted and plotted using GraphPad Prism 8 software.
Results
Overview of the sRNA libraries
To explore and analyse the miRNA expression profiling of S. exigua infected with B. thuringiensis GS57, sRNA libraries were constructed and sequenced for infected (Bt) and uninfected S. exigua larvae (treated with sterilised water) (CK). A total of 608,208,403 raw reads were obtained by using the high-throughput sequencing. After removing low-quality reads, the clean reads ranging from 18 to 30 nt were kept (table S1). The remaining clean reads were 38,305,791 (Bt 1), 18,274,117 (Bt 2), 17,883,388 (Bt 3), 30,017,639 (CK 1), 27,137,717 (CK 2) and 23,193,109 (CK 3), respectively (table 1).
Table 1. sRNA libraries of Bt and CK groups
rRNA, ribosomal RNA; scRNA, small RNA in cytoplasm; snRNA, small nuclear RNA; snoRNA, nucleolar small RNA; tRNA, transport RNA; Repbase, repetitive reads; Unannotated, unannotated reads.
A total of 233 known miRNAs and 726 novel miRNAs were identified in both libraries (table S3). The sequence of identified mature miRNAs and pre-miRNAs were listed in table S4. The lengths of both known and novel miRNAs showed a peak at 22 nt (fig. 1a, b). The percentage of the first base bias towards uracil (U) was 47.19% in known miRNAs and 44.16% in novel miRNAs, respectively (fig. 1c, d), among which 364 miRNAs were divided into 108 families according to the sequence conservation (table S4).
Figure 1. Length distribution and nucleotide bias of known and novel miRNAs. (a) Length distribution of known miRNAs, (b) length distribution of novel miRNAs, (c) nucleotide bias of known miRNAs, (d) nucleotide bias of novel miRNAs.
miRNAs expression and DE miRNAs profiling
The TPM values of identified miRNAs were shown (table S4). The TPM values of pca-bantam-3p, bmo-miR-276-3p and sfr-miR-2766-3p were the highest in Bt and CK RNA libraries (table S4). DE miRNAs were identified and analysed to understand the function of miRNAs in S. exigua infected by B. thuringiensis GS57 strain. As a result, miRNAs with highly similar sequences were identified as the same cluster. Subsequently, 16 up-regulated miRNAs and 34 down-regulated miRNAs were shown in the two RNA libraries (fig. 2a). The highest expressions of sfr-miR-277-3p, pxy-miR-277, mse-miR-277 and hme-miR-277 from CK RNA library were observed; however, the expression levels of novel_miR-514, novel_miR-559, novel_miR-471 and novel_miR-485 identified from Bt RNA library were the highest among DE miRNAs, respectively (fig. 2a and table S4). The expression of novel_miR-571 and miR-123 from Bt RNA library was down-regulated 4.32 folds and up-regulated 5.07 folds compared with the CK RNA library (table S5).
Figure 2. Cluster analysis diagram and verification of the DE miRNAs in Bt and CK groups. (a) Cluster analysis diagram and verification of the DE miRNAs in Bt and CK groups, (b) relative expression of bantam-3p (including mse-bantan, bmo-bantam-3p and hme-bantam), (c) relative expression of miR-277-3p (including sfr-miR-277-3p, pxy-miR-277, hme-miR-277 and mse-miR-277), (d) relative expression of miR-993-5p (including mse-miR-993, bmo-miR-993a-5p and dqu-miR-993-5p), (e) relative expression of novel_miR-170, (f) relative expression of novel_miR-123, (g) relative expression of miR-929-3p (bmo-miR-929-3p). Clustering was performed with log10 (TPM + 1×10−6) values. Columns represent different samples; rows represent different miRNAs. Red blocks represent the higher expressed miRNAs; blue blocks represent the lower expressed miRNAs. The data are shown as the mean ± SE. The differences of DE miRNA expression level between Bt and CK groups were marked with ‘*’ (0.01 < P < 0.05) or ‘**’ (P < 0.01) based on Student's t test.
To verify the accuracy of sRNA sequencing, six DE miRNAs (4 known and 2 novel, 4 up-regulation and 2 down-regulation) (tables S4 and S5) were randomly selected and confirmed the expression level by using RT-qPCR. Analysis of gene relative expression levels showed that bantam-3p (including mse-bantan, bmo-bantam-3p and hme-bantam) (P = 0.0158), miR-277-3p (including sfr-miR-277-3p, pxy-miR-277; hme-miR-277 and mse-miR-277) (P = 0.0316), novel_miR-123 (P = 0.0379) and bmo-miR-929-3p (P = 0.0016) were up-regulated (fig. 2b, c, f and g), whereas miR-939-5p (including mse-miR-993; bmo-miR-993a-5p and dqu-miR-993-5p) (P = 0.0024) and novel_miR-170 (P = 0.0141) were down-regulated in S. exigua infected by B. thuringiensis GS57 (fig. 2d, e).
GO enrichment, KEGG pathway analysis and expression profiling of target genes of DE miRNAs
Total transcripts were used to predict novel potential target genes of DE miRNAs, which was accounted for more potential functions of DE miRNAs. Using GO annotation enrichment, target genes were classified into cellular components, molecular functions and biological processes (fig. 3). KEGG pathway analysis showed that the most target genes are functional in endocytosis and autophagy. Moreover, the remarkable 12 target genes are related to MAPK signalling pathway, five genes are related to ABC transporters and five genes are related to fatty acid metabolism (fig. 4).
Figure 3. Gene ontology (GO) annotation of the target genes of DE miRNAs in Bt and CK groups. x-axis: the GO annotation; left-y-axis: the percentage of genes; right-y-axis: the number of genes.
Figure 4. The most enriched KEGG pathways based on the target genes of DE miRNAs in Bt and CK groups. X-axis: the percentage of annotated genes match to the pathway; y-axis: the pathway names. The different colour of column indicates different types of KEGG pathway.
To understand what kind of biological process do DE miRNAs involve in, the expression level of potential target genes of DE miRNAs were conducted. The selected potential target genes mainly involved in insect development, reproduction and immunity. The results showed that the expression of HF086_004439 (ATP-binding cassette sub-family A member 5-like, ABCA5, GenBank: KAH9638909.1) (P = 0.0004), HF086_017302 (sushi, von Willebrand factor type A, SVWC, GenBank: KAH9638210.1) (P = 0.0141) and HF086_016507 (ATP-binding cassette sub-family D member-like, ABCD, GenBank: KAH9643957.1) (P = 0.0021) were down-regulated after S. exigua was infected by B. thuringiensis GS57. The expression of HF086_003217 (phytanoyl-CoA dioxygenase, peroxisomal-like, PHYD, GenBank: KAH9635463.1) (P = 0.0240) and HF086_006982 (peptidoglycan-recognition protein SB2-like, PGRP, GenBank: KAH9637338.1) (P = 0.0366) were up-regulated after 72 hpi. The expression levels of HF086_011938 (serine/threonine-protein kinase MARK2-like, SPK, GenBank: KAH9632477.1) (P = 0.4755) were almost similar in Bt and CK (fig. 5).
Figure 5. Relative expression of target genes of DE miRNAs between Bt and CK groups. The data are shown as the mean ± SE. The differences of expression level between Bt and CK groups were marked with ‘*’ (0.01 < P < 0.05), ‘**’ (0.001 < P < 0.01), ‘***’ (P < 0.001) or ‘ns’ (no significant difference) based on Student's t test.
Characterisation of miRNA core genes of S. exigua
To identify the miRNA core genes of S. exigua, we obtained amino acid sequences of SeDicer1 and SeAgo1 by BLAST. Conserved domain analysis showed that SeDicer1 protein has one helicase superfamily C-terminal (HELICc) domain, one PiWi-Argonaute-Zwille (PAZ) domain, two ribonuclease III C-terminal (RIBOc) domains and one dsRNA binding motif (DSRM), which is same as S. litura and D. melanogaster (fig. 6a). The SeAgo1 protein has one N-terminal domain of argonaute (ArgoN), one PAZ domain, one Piwi domain and one DSRM domain, which is same as the insects mentioned in fig. 6b. Phylogenetic trees were constructed to examine homologues of SeDicer1 and SeAgo1 in different insects. The SeDicer1 and SeAgo1 were placed with S. litura Dicer1 and Ago1, suggesting a high homologous relationship between S. exigua and S. litura (fig. 6). There were higher homologues in Lepidoptera (fig. 6).
Figure 6. Domain organisation and phylogenetic analysis of Dicer and Ago proteins. (a) Schematic depiction of domain organisation and evolutionary relationships of Dicer proteins, (b) schematic depiction of domain organisation and evolutionary relationships of Ago proteins. HELICc, helicase superfamily C-terminal domain; PAZ, PiWi-Argonaute-Zwille domain; RIBOc, ribonuclease III C-terminal domains; dsRM, double-stranded RNA binding motifs; ArgoN, N-terminal domain of argonaute. For Dicer1, Spodoptera exigua, CAH0698219.1; Drosophila melanogaster, NP_524453.1; Spodoptera litura, XP_022832341.1; Manduca sexta, XP_037296641.1; Tribolium castaneum, XP_008199045.1; Amyelois transitella, XP_013188945.1; Papilio xuthus, KPJ05873.1 were used to analyse homology of Dicer1 proteins in different insects. For Ago1, Spodoptera exigua, KAF9421443.1; Bombyx mori, NP_001095931.1; Drosophila melanogaster, NP_001246314.1; Samia ricini, AID68365.1; Mayetiola destructor, AFX89034.1; Nilaparvata lugens, AGH30326.1; Spodoptera litura, AHC98009.1; Blattella germanica, CCV01212.1; Tribolium castaneum, KYB26000.1 were used to analyse homology of Ago1 proteins in different insects.
Effect of B. thuringiensis GS57 infection on the expression of SeDicer1 and SeAgo1
To verify the effect of B. thuringiensis GS57 infection on the expression of SeDicer1 and SeAgo1, S. exigua was reared on the diet covered with B. thuringiensis GS57. In 24, 48, 72 and 96 hpi, the relative expression level of SeDicer1 and SeAgo1 was determined using RT-qPCR. As shown in fig. 7a, the SeDicer1 expression of S. exigua treated with B. thuringiensis GS57 were higher than those of CK at 72 hpi (P = 0.0002). The relative expression levels of SeDicer1 (P < 0.0001) and SeAgo1 (P = 0.0043) were significantly decreased at 96 h following B. thuringiensis infection comparing with those in the CK (fig. 7a, b). The results suggested that the infection of B. thuringiensis GS57 can affect the expression level of SeDicer1 and SeAgo1 in S. exigua.
Figure 7. Relative expression of SeDicer1, SeAgo1 and effect of dsSeDicer1 on mRNA, miRNAs expression and survival rate of S. exigua infected by B. thuringiensis GS57. (a) Relative expression level of SeDicer1, (b) relative expression level of SeAgo1, (c) relative expression level of SeDicer1 after injecting dsSeDicer1, (d) relative expression levels of miRNAs. (e) Survival rate of S. exigua. The data are shown as the mean ± SE. The data are shown as the mean ± SE. The differences of expression level between different groups were marked with ‘*’ (0.01 < P < 0.05), ‘**’ (0.001 < P < 0.01), ‘***’ (P < 0.001) or ‘ns’ (no significant difference) based on Student's t test. The differences of survival rate of S. exigua between different groups were marked with ‘*’ (P < 0.05) based on Logrank (Mantel-cox test).
Effect of dsRNA on SeDicer1 and miRNAs expression
After injection of dsRNA, the expression of SeDicer1 was reduced by 52% (P = 0.0006), comparing with dsGFP group (fig. 7c). To identify whether SeDicer1 can still effectively regulate miRNA expression when S. exigua was infected by B. thuringiensis GS57, the relative expression of randomly selected DE miRNAs was determined. The relative expression of six DE miRNAs of dsSeDicer1 group was significantly (bantam-3p: P = 0.0004; miR-277-3p: P = 0.0174; miR-993: P = 0.0112; novel_miR-170: P = 0.0118; novel_miR-123: P = 0.0198; miR-929-3p: P = 0.0212) lower than those of dsGFP group (fig. 7d), suggesting that miRNAs were down-regulated after SeDicer1 silencing.
Effect of dsRNA on survival rates of S. exigua
Knockdown of SeDicer1 increased mortality in the larvae following B. thuringiensis GS57 infection, compared with the dsGFP (P = 0.0353) (fig. 7e).
Discussion
Bacillus thuringiensis could produce a variety of insecticidal proteins, which have been widely used in the management of pests. Compared with the insecticidal mechanism, the deep understanding of molecular mechanism of a series of detrimental effects caused by B. thuringiensis infection is limited. These effects were mainly relative to development, reproduction and immunity. For instance, the Cry3Aa protein reduced the growth and reproduction rate of S. littoralis (Hussein et al., Reference Hussein, Habuštová and Sehnal2005). As the previous studies described, B. thuringiensis Cry1C, Cry1Ac and Cry1Ca insecticidal proteins reduced the adult lifespans of Heliothis virescens and S. exigua (Grove et al., Reference Grove, Kimble and McCarthy2001; Zhang et al., Reference Zhang, Ma, Wan, Mu and Li2013b). In addition, the B. thuringiensis 46 isolate did not affect the development of M. separate larvae, but indirectly reduced the hatching rate of the offspring. Adult lifespan and fecundity were decreased when M. separate adult was treated with this strain (Yu et al., Reference Yu, Wang, Sun, Zhang, Li, Zhang and Hou2021). A previous study mainly focused on the regulation of miRNAs in the middle-gut immune response of Plutella xylostella infected by B. thuringiensis HD-73 strain (Li et al., Reference Li, Xu, Zheng, Zheng, Shakeel and Jin2019). For S. exigua, whether and what kinds of miRNAs regulate detrimental effects caused by B. thuringiensis infection on development, reproduction and immunity are still in question.
Here, we constructed sRNA libraries of Bt and CK groups (table 1). In the previous study, a total of 127 known miRNAs were identified from different instar S. exigua larvae (Zhang et al., Reference Zhang, Huang, Yin, Lee, Jia, Liu, Yu, Pennerman, Chen and Guo2015). Our study identified more known and novel miRNAs in S. exigua (table S4), which is helpful for studies of S. exigua miRNAs. The length of known miRNAs and novel miRNAs peaked at 22nt (fig. 1a, b), which conforms to the standard size of animal miRNAs. The first nucleotide of known miRNAs and novel miRNAs was biased towards uracil (U) (fig. 1c, d), which was consistent with first base preference of miRNA. In addition, GO enrichment and KEGG pathway analysis revealed that DE miRNAs caused by B. thuringiensis infection might regulate fatty acid metabolism, MAPK signalling pathway, lysosome and phagosome via potential target genes, which is consistent with B. thuringiensis Cry1Ac-resistant P. xylostella strain (Yang et al., Reference Yang, Xu, Lin, Chen, Lin, Song, Bai, You and Xie2021). Our studies also found that DE miRNAs might involve in the insect hormone biosynthesis, Toll and Imd signalling pathway and FoxO signalling pathway (unpublished data). In our study, total transcripts were used to predict target genes that are related to more pathways (figs 3 and 4), improving the deep understanding of miRNA function in insects infected by B. thuringiensis. However, prediction using total transcripts would result in more incorrect targets, compared with prediction according to 3’ untranslated region (3’ UTR). Thus, more direct evidence of the targeting relationship between miRNAs and potential target genes should be further studied using miRNA mimics and dual luciferase reporter assay (Lee et al., Reference Lee, Kim, Hwang, Jeong, Chung, Lee and Lee2008; Jin et al., Reference Jin, Chen, Liu and Zhou2013; Clément et al., Reference Clément, Salone and Rederstorff2015).
miRNAs are a class of small non-coding RNA, which can involve in regulating a series of biological processes, including insect development, reproduction, immunity and host–pathogen interactions (Hipfner et al., Reference Hipfner, Weigmann and Cohen2002; Ling et al., Reference Ling, Kokoza, Zhang, Aksoy and Raikhel2017; Liu et al., Reference Liu, Zhang, Gao, Feng, Han, Zhang, Bai, Han, Hu, Dai, Zhang and Tong2021). In our study, DE miRNAs were caused by B. thuringiensis infection, and might be critical regulatory factors in the interaction between B. thuringiensis and S. exigua. The role of DE miRNAs mainly depended on the function of potential target genes. The function of potential target genes was predicted in development, reproduction and immunity (figs 3 and 4). Based on this, the relative expression of the five genes from the selected six potential target genes, which are related to insect development, reproduction and immunity, showed differences between Bt and CK (fig. 5). Among them, ABCA5 plays an important role in lysosome-related immunity (Kubo et al., Reference Kubo, Sekiya, Ohigashi, Takenaka, Tamura, Nada, Nishi, Yamamoto and Yamaguchi2005), and PGRP can recognise pathogens and activate the immune response. Furthermore, ABCD and PHYD expression might cause peroxisome-related immunity (Petriv et al., Reference Petriv, Pilgrim, Rachubinski and Titorenko2002), implying miRNAs could be involved in immune response to B. thuringiensis via target genes. In addition, ABCD can be involved in the transformation of fatty acids and acyl-CoAs (Theodoulou et al., Reference Theodoulou, Holdsworth and Baker2006; Tian et al., Reference Tian, Song, He, Zeng, Xie, Wu, Wang, Zhou and Zhang2017), suggesting DE miRNAs might regulate the development of S. exigua that were infected by B. thuringiensis GS57. The previous study showed that the expression level of SVWC decreased as B. mori was infected by fungus, and directly regulated moulting and development, both demonstrating that SVWC can be involved in insect development and immunity (Han et al., Reference Han, Lu, Yuan, Huang, Beerntsen, Huang and Ling2017). The expression of SVWC of S. exigua infected by B. thuringiensis GS57 was decreased. Although the function of SVWC in immunity might be caused by the way of fungal infection, it is also noteworthy. More potential target genes were predicted which proved the possibility that DE miRNAs were involved in other biological processes of S. exigua infected by B. thuringiensis GS57. These results suggested that miRNA could regulate multiple processes of insects infected by B. thuringiensis via target genes.
The formation and function of miRNAs mainly depend on two miRNA relative core genes, Dicer1 and Ago1 (Llave et al., Reference Llave, Xie, Kasschau and Carrington2002; Filipowicz et al., Reference Filipowicz, Bhattacharyya and Sonenberg2008). A previous study proved that the expressions of Dicer1 and Ago1 were affected in H. armigera injected with B. thuringiensis (Baradaran et al., Reference Baradaran, Moharramipour, Asgari and Mehrabadi2019), which is consistent with this research. During this period of infection, the expression of Ago1 significantly up-regulated at 24 and 72 hpi (Baradaran et al., Reference Baradaran, Moharramipour, Asgari and Mehrabadi2019), which is different from our study (fig. 7). The difference might be caused by infection method. The infection of B. thuringiensis caused by direct bacteria injection would lead to a sharp increase of genes expression. Here, we found that the expression of SeDicer1 in Bt groups was higher than those of CK groups in 72 hpi, suggesting that B. thuringiensis infection directly affected the expression of miRNA core genes, and then indirectly affected the formation of miRNAs. Thus, B. thuringiensis infection led to the change of expression of DE miRNAs. The up-regulation of selected potential target genes, SePGRP and SePHYD, might suggest that the immune response was activated to defend against B. thuringiensis infection. The expression of SeABCD decreased, which affected the fatty acid metabolism of S. exigua, mitigating the fitness costs associated with the defence against the B. thuringiensis infection. Thus, these biology processes were driven by SeDicer1 and SeAgo1 via miRNAs.
Recent studies demonstrated that the expression of Let-7 and miR-184 decreased after Dicer1 knockdown (Rahimpour et al., Reference Rahimpour, Moharramipour, Asgari and Mehrabadi2019; Bidari et al., Reference Bidari, Fathipour, Asgari and Mehrabadi2022). In our study, the silencing of SeDicer1 could lead to the decrease of multiple miRNAs, and promote mortality following B. thuringiensis GS57 infection (fig. 7d, e). The previous studies demonstrated that silencing of Dicer1 led to larvae development arrested or mortality (Zhang et al., Reference Zhang, Lu and Zhou2013a; Rahimpour et al., Reference Rahimpour, Moharramipour, Asgari and Mehrabadi2019; Chen et al., Reference Chen, Yang, Yang, Liu, Wang, Luo, Tang, Yi, Huang, Liu and Liu2023). For example, the knockdown of Dicer1 decreased the survival rate in Sogatella furcifera (Zeng et al., Reference Zeng, Long, Yang, Zhou, Yang, Wang and Jin2023). These results proved that Dicer1 played a key role in insect growth and development. Thus, silencing of SeDicer1 might disturb the larvae development, and impair the synthesis of miRNAs to weaken the defence ability of S. exigua to B. thuringiensis infection, which both led to the increased mortality of S. exigua. However, the previous study reported that silencing of Dicer1 could decrease the replication of B. thuringiensis bacteria, and indirectly decrease the mortality rate of H. armigera injected by B. thuringiensis (Baradaran et al., Reference Baradaran, Moharramipour, Asgari and Mehrabadi2019), which is different from this study. The B. thuringiensis infection caused by injection is mainly related to septicaemia, but feeding-infection might be related to multiple damages. This difference of infection method might affect the function of Dicer1. In addition, miRNAs should load onto Ago1 protein to regulate the post-transcription of target genes. The effect of silencing of Ago1 on the expression of miRNAs should be still investigated in the future. The previous studies have proved that insects would develop resistance under multi-generation B. thuringiensis stress (Tabashnik et al., Reference Tabashnik, Carrière, Dennehy, Morin, Sisterson, Roush, Shelton and Zhao2003, Reference Tabashnik, Van Rensburg and Carrière2009). Therefore, multiple management strategies should be explored to manage pests and slow the development of B. thuringiensis resistance. RNAi-based pest management is being developed and exploited, which mainly focuses on the expression of exogenous dsRNA in transgenic plants and spray-induced gene silencing insecticide (Price and Gatehouse, Reference Price and Gatehouse2008; Mezzetti et al., Reference Mezzetti, Smagghe, Arpaia, Christiaens, Dietz-Pfeilstetter, Jones, Kostov, Sabbadini, Opsahl-Sorteberg, Ventura, Taning and Sweet2020). For example, the transgenic Nicotiana tabacum with the expression of H. armigera HR3 dsRNA could disrupt H. armigera development (Xiong et al., Reference Xiong, Zeng, Zhang, Xu and Qiu2013). Our results demonstrated that the knockdown of SeDicer1 would enhance the insecticidal toxicity of B. thuringiensis (fig. 7e). Thus, transgenic plants and spray-induced SeDicer1 silencing along with B. thuringiensis application might be a new potential strategy.
Traditional RNAi-based pest management mainly knockdowns specific target genes using dsRNA or siRNA. Recently, more pest management is focused on the miRNA-based RNAi. In a recent study, miR-34-5p has been considered as a novel molecular target for Lepidoptera pests. The over-expression or knockdown of miR-34-5p can lead to high mortality, low fecundity and developmental defects (Li et al., Reference Li, Zhu, Sun, Zheng, Liang and Gao2022a). In addition, the previous study demonstrated that engineering strains of Beauveria bassiana over-expressed exogenous miRNA-8 or miRNA-375 were constructed and had higher fungal efficacy in controlling pests (Cui et al., Reference Cui, Wang, Li, Sun, Jiang, Zhou, Liu and Wang2022). However, the latest study showed that insects could compensate for the knockdown of specific miRNAs by maintaining the expression of target genes by other pathways (Zuo et al., Reference Zuo, Wang, Ren, Pei, Aioub and Hu2022), revealing the screening difficulty of key miRNAs. These findings suggested that miRNA-based pest management can control pests, suggesting a possibility of these novel molecular targets. In the future, the insecticidal potential of DE miRNAs should be identified and explored to control S. exigua.
Conclusion
In this study, we found that B. thuringiensis infection caused changes in miRNAs and miRNA core genes. The KEGG pathway analysis revealed that potential target genes of DE miRNAs were associated with ABC transporters, fatty acid metabolism and MAPK signalling pathways, suggesting miRNAs might be involved in regulating the change of development, reproduction and immunity caused by B. thuringiensis infection. These findings help to reveal the non-traditional effects of B. thuringiensis on insects in terms of miRNAs. In addition, the expression of SeDicer1 could directly regulate the synthesis of miRNA, and affect the insecticidal activity of B. thuringiensis. In conclusion, these results may have important implications for reconsidering the functions of miRNAs and miRNAs core genes in S. exigua infected by B. thuringiensis, and exploring miRNAs-based pest management.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485324000300.
Acknowledgements
This work was supported by the earmarked fund for Modern Agro-industry Technology Research System (CARS-13), the National Natural Science Foundation of China (Grant No. 31471775) and the Foundation of the Graduate School of the Chinese Academy of Agricultural Sciences (CAAS) (1610042022005).
Author contributions
Bo Gao: original manuscript writing, conceptualisation, methodology and investigation. Yu-Jie Ji: methodology and investigation. Dan Zhao: methodology and investigation. Lu Zhang: formal data analysis and software. Han Wu: formal data analysis and software. Yi-Fan Xie: formal data analysis and software. Qiu-Yu Shi: formal data analysis and software. Wei Guo: review and editing, conceptualisation, funding acquisition.
Competing interests
None.
