Determination of internal controls for quantitative gene expression of Spodoptera litura under microbial pesticide stress

TL;DR

This study identifies the most stable reference genes for gene expression analysis in Spodoptera litura under different microbial pesticide treatments.

Contribution

The study provides experimentally validated reference genes for qRT-PCR normalization in S. litura under pesticide stress.

Findings

RPLP0 and ACT5C are the best reference genes for direct pesticide treatment.

RPS13 and GAPDH are most stable for comprehensive treatment conditions.

RPS13 and RPLP0 are the most reliable across all treatment types.

Abstract

Quantitative real-time polymerase chain reaction (qRT-PCR) has become a commonly used method for the quantification of gene expression. However, accurate qRT-PCR analysis requires a valid internal reference for data normalization. To determine the valid reference characterized with low expression variability among Spodoptera litura samples after microbial pesticide treatments, nine housekeeping genes, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), arginine kinase, ubiquitin C, actin-5C (ACT5C), actin, ribosomal protein S13 (RPS13), tubulin, acidic ribosomal protein P0 (RPLP0) and ubiquinol-cytochrome c reductase, were evaluated for their suitability using geNorm, Normfinder, BestKeeper, RefFinder and the comparative delta CT methods in this study. S. litura larvae after direct treatment (larvae were immersed in biopesticides), indirect treatment (larvae were fed with biopesticide immersed artificial diets) and comprehensive treatment (larvae were treated with the first two treatments in sequence), respectively with Metarhizium anisopliae, Empedobacter brevis and Bacillus thuringiensis, were investigated. The results indicated that the best sets of internal references were as follows: RPLP0 and ACT5C for direct treatment conditions; RPLP0 and RPS13 for indirect treatment conditions; RPS13 and GAPDH for comprehensive treatment conditions; RPS13 and RPLP0 for all the samples. These results provide valuable bases for further genetic researches in S. litura.

Figures

Regression lines for the nine candidate reference genes. Each dot indicates the mean of duplicate samples (n = 3). The most stable reference gene has the closest fit to the regression line.Figure 2Scatterplot of residuals analysis. The residuals of CT values were evaluated by linear regression in Fig. [1](#Fig1) and the difference between the actual value and the calculated value for each gene.

The average expression stability value M and the pairwise variation V of candidate genes as determined by geNorm analysis. (a and b) direct treatment conditions; (c and d) indirect treatment conditions; (e and f) comprehensive treatment conditions; (g and h) all treatments.Figure 4Expression stabilities of the candidate reference genes as determined by Normfinder software. a direct treatment conditions; b indirect treatment conditions; c comprehensive treatment conditions; d all treatments.Table 1Stability of candidate reference genes as determined by BestKeeper software.GeneDirect treatment c

The relative expression levels of C-Type Lectin genes in *S. litura* after 24 h and 72 h of treatment with microbial pesticide. (a and d) direct treatment conditions; (b and e) indirect treatment conditions; (c and f) comprehensive treatment conditions; (a–c) *SlCTL* expression levels when normalized by the internal references determined from different treatment conditions; (d-f) *SlCTL* expression levels when normalized by the internal references determined from all treatments. Values are expressed as mean ± SE (n = 3). Different letters above each bar indicate statistical differences, as det

The relative expression levels of C-Type Lectin genes in *S. litura* after 24 h and 72 h of treatment with microbial pesticide. (a and d) direct treatment conditions; (b and e) indirect treatment conditions; (c and f) comprehensive treatment conditions; (a–c) *SlCTL* expression levels when normalized by the internal references determined from different treatment conditions; (d-f) *SlCTL* expression levels when normalized by the internal references determined from all treatments. Values are expressed as mean ± SE (n = 3). Different letters above each bar indicate statistical differences, as det

Full text

Introduction

Spodoptera litura Fabricius (Lepidoptera: Noctuidae) is a globally distributed polyphagous pest that damages approximately 389 species of plants including vegetables, fruits, cotton, and tobacco; the most commonly affected plants are crop species^1^. The high reproductive potential of this species and intense nutritional requirements of its larva means that most damage is incurred over a short period of time^2–6^. Outbreaks of S. litura have been reported in many countries, including China, India, Pakistan, Japan, Indonesia and Australia^7–13^. Over the past few decades, chemical control has been utilized as the main strategy for managing S. litura. However, the development of resistance to chemical pesticides in this species leads to subsequent management failure, posing a serious threat to global agricultural production^14–17^.

Over recent years, more biological approaches have been developed to effectively control S. litura^18,19^. Biopesticides are the biological agents that are used to control pests, and are derived from fungi, bacteria, viruses, plants, animals, and certain minerals. Of all biopesticides, microbial pesticides are becoming increasingly more important and have gained significant popularity because they are safe and environmentally friendly^20,21^. The most widely used microbial pesticides are strains of Bacillus, Metarhizium and Empedobacter, these have all been demonstrated to be effective against various Lepidoptera pests^22–28^. However, in a manner similar to that of chemical control methods, some target herbivores, including S. litura, have developed resistance to microbial pesticides, including B. thuringiensis^29–33^.

Molecular technologies, especially quantitative real-time polymerase chain reaction (qRT-PCR), have been used extensively in genetic studies relating to the mechanisms of immunity in insects^34–36^. Significant changes in the levels of gene expression can reflect biological changes in insects under different experimental conditions. To investigate the specific changes in immune-associated genes in S. litura under different microbial pesticide stress conditions, especially under different exposure treatments to pesticide, valid internal references for qRT-PCR analyses were screened. One group of S. litura larvae was treated directly with M. anisopliae, E. brevis and B. thuringiensis, respectively, by immersion in biopesticide; another group of larvae was indirectly treated by feeding the larvae with artificial diets immersed in biopesticide; the final group of larvae was treated using the first two treatment modes in sequence; we referred to this strategy as the comprehensive treatment. Nine housekeeping genes, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), arginine kinase (AK), ubiquitin C (UBC), actin-5C (ACT5C), actin (ACT), ribosomal protein S13 (RPS13), tubulin (TUB), acidic ribosomal protein P0 (RPLP0) and ubiquinol-cytochrome c reductase (UCCR), in samples of S. litura following different treatments with microbial pesticide, were evaluated for their suitability as normalization references using geNorm, Normfinder, BestKeeper, RefFinder and the comparative delta CT methods.

Results

Qualities of total RNA

Total RNAs were extracted from S. litura after 6, 12, 24, 48 and 72 h of direct treatments, and after 24, 48, 72 h of indirect treatments and comprehensive treatments, respectively. The concentrations and purities of the total RNA isolated from S. litura samples were determined with a GeneQuant Pro RNA/DNA Calculator (GE Healthcare, Piscataway, NJ, USA). The total RNA concentrations ranged from 711.3 to 1654.8 ng·μL^−1^ for the directly treated samples, from 810.9 to 1674.3 ng·μL^−1^ for the indirectly treated samples, and from 767.9 to 1284.1 ng·μL^−1^ for the comprehensive groups. The A260/A280 ratios ranged from 2.04 to 2.11 for the directly treated samples, from 1.91 to 2.11 for the indirectly treated samples, and from 1.93 to 2.12 for the comprehensive groups (Supplementary Table S1). The integrity of all total RNA samples was confirmed by 1.0% agarose gel electrophoresis.

PCR amplification efficiencies

For each primer pair, the single peak melting curves indicated that a unique product was amplified (Supplementary Fig. S1). The products were sequenced and BLAST searches were performed at http://www.ncbi.nlm.nih.gov/blast/. BLASTn revealed that the products had 100% identity with the fragment sequences on which the primer design as based. The PCR amplification efficiency and the coefficient of determination (R^2^) were 99.40 and 0.9971 for GAPDH, 109.39 and 0.9972 for AK, 109.53 and 0.9998 for UBC, 94.00 and 0.9994 for ACT5C, 96.86 and 0.9984 for ACT, 106.41 and 0.9977 for RPS13, 91.22 and 0.9985 for TUB, 102.18 and 0.9968 for RPLP0, 106.64 and 0.9998 for UCCR, respectively (Supplementary Fig. S2).

Expression profiles of the candidate reference genes

According to the results of crude expression levels and stability of each gene from our previous research on S. litura transcriptome under microbial pesticide stress, the nine housekeeping genes (GAPDH, AK, UBC, ACT5C, ACT, RPS13, TUB, RPLP0 and UCCR) were chosen to serve as the candidate reference genes for this study. The expression levels of the nine candidate reference genes in S. litura samples were investigated using a SYBR Green-based qPCR assay which was performed in triplicate. The entire experiment was then repeated. The mean cycle threshold (CT) values ranged from 15.35 (ACT) to 26.50 (UBC) in all samples, from 15.37 to 20.48 for ACT5C, from 15.35 to 22.41 for ACT, from 18.43 to 22.83 for GAPDH, from 19.30 to 22.76 for RPLP0, from 17.92 to 22.46 for RPS13, from 18.96 to 24.56 for TUB, from 22.99 to 26.50 for UBC, from 21.80 to 26.38 for UCCR, and from 17.23 to 21.95 for AK (Fig. 1). The CT values in highest to lowest order were as follows: UBC, UCCR, TUB, RPLP0, RPS13, GAPDH, AK, ACT and ACT5C. The residuals of CT values were evaluated by linear regression and the difference between the actual value and the calculated value for each gene (Fig. 1); CT values were ranked as follows (highest to lowest stability): RPLP0, UBC, GAPDH, RPS13, UCCR, AK, TUB, ACT5C and ACT, according to the distributions of residuals (Fig. 2).Figure 1. Regression lines for the nine candidate reference genes. Each dot indicates the mean of duplicate samples (n = 3). The most stable reference gene has the closest fit to the regression line.Figure 2. Scatterplot of residuals analysis. The residuals of CT values were evaluated by linear regression in Fig. 1 and the difference between the actual value and the calculated value for each gene.

Analysis of gene expression stability

Direct treatment

According to geNorm analysis, the stability rankings from the most stable to the least stable gene were ACT, ACT5C, RPLP0, UBC, TUB, RPS13, UCCR, GAPDH and AK (Fig. 3a). Furthermore, geNorm analysis revealed that the pair-wise variation value V5/6 was below the proposed 0.15 cut-off (Fig. 3b). This result suggested that the average of the top five genes would be the optimal normalization factor for further experiments. Normfinder analysis identified RPLP0, ACT5C and UBC as the most stable genes (Fig. 4a). According to the standard deviation (SD) and coefficient of variation (CV) values in Table 1, BestKeeper analysis identified RPLP0, UBC and UCCR as the most stable genes. The stability rankings generated by the delta CT and RefFinder methods identified RPLP0, ACT5C and UCCR as the most stable genes (Fig. 5a,e). Moreover, all software tools identified AK as the least stable gene while most software tools identified RPLP0 and ACT5C as the top two most stable genes.Figure 3. The average expression stability value M and the pairwise variation V of candidate genes as determined by geNorm analysis. (a and b) direct treatment conditions; (c and d) indirect treatment conditions; (e and f) comprehensive treatment conditions; (g and h) all treatments.Figure 4. Expression stabilities of the candidate reference genes as determined by Normfinder software. a direct treatment conditions; b indirect treatment conditions; c comprehensive treatment conditions; d all treatments.Table 1. Stability of candidate reference genes as determined by BestKeeper software.GeneDirect treatment conditionsIndirect treatment conditionsComprehensive treatment conditionsAll samplesRankSDCVrRankSDCVrRankSDCVrRankSDCVr**ACT5C50.744.320.86860.633.690.89881.085.940.94670.885.020.910ACT40.713.990.77280.724.110.87591.357.220.93460.864.790.857GAPDH80.944.690.75840.582.900.65830.622.960.93540.834.110.813RPLP010.512.510.87420.502.380.81220.542.530.74320.622.990.842RPS1360.804.010.84610.422.090.74040.643.030.94230.753.720.883TUB70.864.200.82830.522.510.60450.783.430.85191.125.280.810UBC20.532.150.71650.612.440.62710.532.120.31710.572.320.550UCCR30.702.960.81690.793.100.71670.913.760.94580.953.910.617AK90.984.980.80970.643.310.69060.904.510.84450.844.300.766Figure 5Expression stabilities of the candidate reference genes as determined by the comparative delta CT method (a–d) and RefFinder software (e–h). (a and e) direct treatment conditions; (b and f) indirect treatment conditions; (c and g) comprehensive treatment conditions; (d and h) all treatments.

Indirect treatment

According to geNorm analysis, the stability rankings from the most stable to the least stable gene were RPS13, RPLP0, TUB, GAPDH, ACT, ACT5C, AK, UBC and UCCR (Fig. 3c). Furthermore, geNorm analysis revealed that the pair-wise variation value V5/6 was below the proposed 0.15 cut-off (Fig. 3d). This result suggested that the average of the top five genes would be the optimal normalization factor for further experiments. Normfinder and delta CT methods identified RPLP0, RPS13 and ACT as the most stable genes (Figs. 4b, 5b). According to the standard deviation and CV values in Table 1, BestKeeper software identified RPS13, RPLP0 and TUB as the most stable genes; these were the same as those identified by geNorm. RefFinder software identified RPLP0, RPS13 and TUB as the top three most stable genes (Fig. 5f); these results were similar to those derived from geNorm and BestKeeper. Moreover, all of the software tools identified UCCR as the least stable gene, and identified RPLP0 and RPS13 as the top two most stable genes.

Comprehensive treatment

According to geNorm analysis, the stability rankings from the most stable to the least stable gene were RPS13, GAPDH, UCCR, RPLP0, TUB, AK, ACT5C, ACT and UBC (Fig. 3e). Furthermore, geNorm analysis revealed that the pair-wise variation value V4/5 was below the proposed 0.15 cut-off (Fig. 3f). This result suggested that the average of the top four genes would be the optimal normalization factor for further experiments. Normfinder and delta CT methods also identified RPS13, GAPDH and UCCR as the top three most stable genes (Figs. 4c, 5c); these results were the same as those generated by geNorm. According to the standard deviation and CV values in Table 1, BestKeeper software identified UBC, RPLP0 and GAPDH as the most stable genes. It’s notable that UBC was identified as the least stable gene by geNorm. RefFinder analysis identified RPS13, GAPDH and RPLP0 as the most stable genes (Fig. 5g). Moreover, all software tools, except for geNorm, identified ACT as the least stable gene. Most of the programs identified RPS13 and GAPDH as the top two most stable genes.

All treatments

Next, we investigated the stability rankings of the nine candidate reference genes for all treated S. litura samples. According to geNorm analysis, the stability rankings from the most stable to the least stable gene were RPS13, RPLP0, ACT5C, ACT, GAPDH, AK, UBC, TUB and UCCR (Fig. 3g). In addition, geNorm analysis revealed that the pair-wise variation value V6/7 was below the proposed 0.15 cut-off (Fig. 3h). This result indicated that the average of the top six genes would be the optimal normalization factor for further experiments. All software packages, except for BestKeeper, identified RPS13, RPLP0 and ACT5C as the most stable genes (Figs. 4d, 5d,h). UBC, RPLP0 and RPS13 were the most stable genes generated by the BestKeeper software (Table 1). Moreover, all software packages, except for BestKeeper, identified UCCR as the least stable gene.

Validation of the candidate internal reference genes

Previous research has shown that C-type lectins (CTLs) participate in pathogen recognition in insects and play diverse roles in a range of immune responses, including opsonization, agglutination, nodule formation, encapsulation, phagocytosis, melanization, prophenoloxidase activation and homeostatic maintenance of the gut microbiome^36–38^. In addition, the expression levels of CTLs in S. litura were demonstrated to undergo change following fungal or bacterial infections^36^. In the present study, we evaluated the relative expression levels of SlCTL in S. litura after 24 h and 72 h of microbial pesticide treatments, using the multiple genes recommended by geNorm, the top three and the top two most stable genes, and the least stable gene identified by most of software packages, as internal references for data normalization. Figure 6 shows that there was no significant difference among the relative expression levels of SlCTL normalized by the following combinations: ACT + ACT5C + RPLP0 + UBC + TUB, RPLP0 + ACT5C + UBC, RPLP0 + ACT5C + UCCR and RPLP0 + ACT5C for the direct treatment (P > 0.05) (Fig. 6a); RPS13 + RPLP0 + TUB + GAPDH + ACT, RPLP0 + RPS13 + ACT, RPLP0 + RPS13 + TUB and RPLP0 + RPS13 for the indirect treatment (P > 0.05) (Fig. 6b) ; RPS13 + GAPDH + UCCR + RPLP0, RPS13 + GAPDH + UCCR, RPS13 + GAPDH + RPLP0 and RPS13 + GAPDH for the comprehensive treatment (P > 0.05) (Fig. 6c). Moreover, the expression levels of SlCTL when normalized by the least stable genes (AK, UCCR, ACT and UBC) under different treatments were significantly different from the levels normalized by the other gene combinations (P < 0.05) (Fig. 6a–c). When take all the treatment samples into consideration, the gene combinations, RPS13 + RPLP0 + ACT5C + ACT + GAPDH + AK, RPS13 + RPLP0 + ACT5C and RPS13 + RPLP0, as well as the least stable gene UCCR, were used for normalizing the expression levels of SlCTL. Figure 6d,f show that the relative expression levels of SlCTL normalized by UCCR were significantly different from the levels normalized by the other gene combinations (P < 0.05). These results demonstrated that using the top two most stable genes was valid for the normalization of SlCTL expression levels under the experimental conditions used herein, and that the use of inappropriate internal references would lead to inaccurate experimental results.Figure 6. The relative expression levels of C-Type Lectin genes in S. litura after 24 h and 72 h of treatment with microbial pesticide. (a and d) direct treatment conditions; (b and e) indirect treatment conditions; (c and f) comprehensive treatment conditions; (a–c) SlCTL expression levels when normalized by the internal references determined from different treatment conditions; (d-f) SlCTL expression levels when normalized by the internal references determined from all treatments. Values are expressed as mean ± SE (n = 3). Different letters above each bar indicate statistical differences, as determined by ANOVA followed by Duncan’s multiple range test (P < 0.05).

Discussion

Housekeeping genes, such as 18S ribosomal RNA (18S rRNA), elongation factor 1 alpha (EF1-a), ribosomal protein L18 (RPL18), ribosomal protein S18 (RPS18), beta actin (β-actin), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), alpha tubulin (α-Tub), beta tubulin (β-Tub), TATA-Box binding protein (TBP), and glucose-6-phosphate dehydrogenase (G6PDH), are commonly used as internal reference genes in insect studies. Of these, ACT, RPS18 and GAPDH have previously been identified as the most stable genes in Apis mellifera following infection with Escherichia coli^39^. RPS3, RPS18, and RPL13a have been identified as the most stable genes in Tribolium castaneum following infection with Beauveria bassiana^40^. For gene studies involving S. litura, the β-actin gene was previously used as an internal reference to determine target gene expression patterns under zinc and Nomuraea rileyi stress^34,41^. EF-1 has been used to normalize expression levels following infection with N. rileyi, SpltNPV and B. thuringiensis, respectively^36^, and other reference genes were identified across different biotic and abiotic experimental conditions^35,42^.

Previous studies on S. litura indicated that some commonly used housekeeping genes exhibited significant variation in expression under different experimental conditions^35,42^. To investigate the molecular mechanisms of immunity in S. litura under different microbial pesticide stress conditions, especially in terms of different treatment modes of exposure to pesticide, we validated a range of internal reference genes for data normalization in the present study. Based on our previous research on the S. litura transcriptome under microbial pesticide stress, the nine housekeeping genes, GAPDH, AK, UBC, ACT5C, ACT, RPS13, TUB, RPLP0 and UCCR, with appropriate FPKM values (fragments per kilo base of transcript per million fragments mapped), were selected to serve as candidate reference genes. Since an appropriate expression level (a CT value between 15 and 30) is important for analyzing internal reference genes^43–47^, the CT values of the nine candidate reference genes were confirmed in all of the treated S. litura samples, and the results showed that the average CT values ranged from 15.35 (ACT) to 26.50 (UBC) (Fig. 1).

The stabilities of the nine candidate reference genes were analyzed by geNorm, Normfinder, BestKeeper, RefFinder and the comparative delta CT methods. The results indicated that in most cases the top three most stable genes ranked by BestKeeper under each treatment condition were different from those obtained from other software packages. The relative expression levels of SlCTL, when normalized by UBC under comprehensive treatment conditions, were significantly different from the levels normalized by the other references (P < 0.05) (Fig. 6c), in particular, UBC was identified as the most stable gene by BestKeeper software (Table 1), but was identified as the least stable gene when ranked by geNorm (Fig. 3e). Consequently, the internal reference genes recommended by BestKeeper should not be applied for further analyses under the experimental conditions described herein. Furthermore, the geNorm manual states that the application of the three best reference genes is a valid normalization strategy in most cases. Therefore, the relative expression levels of SlCTL were also normalized by the combination of the top three most stable genes under each treatment condition. Figure 6 showed that all of the detected gene combinations were valid for the normalization of SlCTL expression levels under the experimental conditions used herein; there was no significant difference between the SlCTL expression levels normalized by the combinations of the top three and the top two most stable genes (P > 0.05). Finally, the preferable reference genes across different treatment conditions according to our overall analysis were as follows: RPLP0 and ACT5C for direct treatment conditions; RPLP0 and RPS13 for indirect treatment conditions; RPS13 and GAPDH for comprehensive treatment conditions; along with RPS13 and RPLP0 for all samples.

RPLP0 is located in the 60S ribosomal subunit, which plays a role in the association of elongation factors with the ribosome during protein synthesis^48–50^, DNA repair^51^, gene expression regulation^52^ and O_2_ consumption cycles^53^. RPS13 is located in the 40S ribosomal subunit and plays a role in peptide chain elongation and translocation of the mRNA:tRNA complex^54^. These two ribosomal proteins were alternately ranked as the most stable gene in S. litura under the experimental conditions used herein. As with studies performed on other insect species, ribosomal protein genes have always been validated as internal controls for qRT-PCR^39,42,55–61^.

In summary, the internal controls applied for qRT-PCR studies on S. litura under different microbial pesticide stress conditions, especially with different treatment modes of exposure to pesticide, were different. The selection of a valid normalization gene is particularly important for experimental reliability. Our findings provide valuable bases for further research on genes in S. litura.

Materials and methods

Insects and biopesticides

S. litura were collected as larvae from cabbage fields at the National Center for Vegetable Improvement (Chongqing, China) in October 2022. The larvae were identified as S. litura by analyzing larval morphology, especially according to the two subtriangular dark spots on each segment of the larva, except for the prothorax. Ten generations of S. litura were reared in the laboratory to minimize the effects of field environments. Larvae were reared on a soybean and wheat bran-based artificial diet^62^. Insects were kept at 26 °C, with a 12 h photoperiod, and a relative humidity of 70%. Three biological pesticide products, 8.0 × 10^9^ spores·mL^-1^ of Metarhizium anisopliae CQMa421 OD (Chongqing Julixin Bioengineering Co., Ltd, Chongqing, China), 1.0 × 10^10^ spores·mL^-1^ of Empedobacter brevis GXW15-4 SC (Zhenjiang Runyu Biological Science and Technology Development Co., Ltd, Jiangsu, China), and 16,000 IU·mg^-1^ of Bacillus thuringiensis WP (King Biotec Corp. Hubei, China), were used in this study. Since our preliminary experiments showed that 50% of S. litura larvae could survive for 7 days after immersion treatment with 2.00 × 10^7^ spores·mL^-1^ of M. anisopliae, 6.25 × 10^7^ spores·mL^-1^ of E. brevis and 2.50 mg·mL^-1^ of B. thuringiensis, respectively, these concentrations of biopesticides were used for subsequent experiments to ensure sufficient time and larvae before they pupated.

Biopesticide stress

Direct treatment

Ten fourth instar day 1 larvae were immersed in M. anisopliae (2.00 × 10^7^ spores·mL^-1^), E. brevis (6.25 × 10^7^ spores·mL^-1^) and B. thuringiensis (2.50 mg·mL^-1^), respectively, for 10 s and allowed to air dry at room temperature. The larvae immersed in distilled water were used as control. Three replicates were prepared for each sample. The larvae were collected at 6, 12, 24, 48 and 72 h post-infection.

Indirect treatment

Ten fourth instar day 1 larvae were fed with 1.0 cm^3^ of an artificial diet immersed in M. anisopliae (2.00 × 10^7^ spores·mL^-1^), E. brevis (6.25 × 10^7^ spores·mL^-1^) and B. thuringiensis (2.50 mg·mL^-1^), respectively, and the larvae fed with distilled water immersed artificial diets were used as control. Untreated artificial diets were provided when the immersed diets were exhausted. Three replicates were prepared for each sample. Since it would take around 24 h for the larvae to finish the biopesticides-immersed artificial diets, the larvae after 24, 48 and 72 h treatment were collected, respectively.

Comprehensive treatment

Ten fourth instar day 1 larvae were immersed in M. anisopliae, E. brevis and B. thuringiensis, and fed with 1.0 cm^3^ of an artificial diet immersed in biopesticide, respectively. The larvae and diets treated with distilled water were used as control. Untreated artificial diets were provided when the immersed diets were exhausted. Three replicates were prepared for each sample. Since it would take around 24 h for the larvae to finish the biopesticides-immersed artificial diets, the larvae after 24, 48 and 72 h treatment were collected, respectively.

Total RNA extraction and cDNA synthesis

Total RNA was isolated from five larvae with TRIzol reagent (Invitrogen, Carlsbad, CA, USA). First-strand cDNAs were then synthesized with PrimeScript™ RT Master Mix (Perfect Real Time) (Takara). The reactions were performed in accordance with the manufacturer’s instructions. The nucleotide sequences of the primers used to amplify cDNA fragments of GAPDH, AK, UBC, ACT5C, ACT, RPS13, TUB, RPLP0 and UCCR for qPCR are shown in Table 2. PCR reactions were conducted in a S1000™ Thermal Cycler (Bio-Rad). The PCR products were then evaluated by 1.0% agarose gel electrophoresis. Bands of the expected sizes were excised and then each fragment was purified using a MiniBEST Agarose Gel DNA Extraction Kit Ver.3.0 (Takara).Table 2. Primer pairs used for qPCR.Gene name (Abbreviation)Accession NoSequence (5′-3′)Product length (bp)PCR efficiencies (%)R^2^Glyceraldehyde-3-phosphate dehydrogenase (GAPDH)HQ012003F: CTGATGCTCCCATGTTCGTG20299.40.9971R: CCAGAGGGTCCATCAACAGTArginine kinase (AK)HQ840714F: GACCTTCTTGGTATGGTGCAATGAAG175109.40.9972R: GTTGGTAGGGCAGAAAGTGAGGAUbiquitin C (UBC)XM_022968576F: CCTTGACGGGTAAAACTATTACGCTTGAA206109.50.9998R: GCCACGGAGACGCAGAACAActin-5C (ACT5C)XM_022975529F: CGAGAAATCGTGCGTGACAT24494.00.9994R: CGTCGCACTTCATGATGGAGActin (ACT)XM_022981497F: CACCTTCTACAACGAGCTGC17496.90.9984R: CCAGAGGCGTACAGAGAGAGRibosomal protein S13 (RPS13)XM_022972107F: GTATGCACGCACCTGGTAAG185106.40.9977R: TCTGACTTGTGCGACTCCATTubulin (TUB)XM_022962426F: GACAACGAGGCCCTATACGA20591.20.9985R: CGAATCCGGGCATGAAGAAGAcidic ribosomal protein P0 (RPLP0)XM_022968776F: GGAAACCAACCCAGCTCTTG198102.20.9968R: GTCTTCTCAGGACCCAGACCUbiquinol-cytochrome c reductase (UCCR)XM_022974041F: GGGCAATTCTCTTTTCATCTCACCCA227106.60.9997R: CACCCATTTCTTTCCTCAAATCTCCACC

Quantitative real-time PCR

qRT-PCR was performed using a qTOWER^3^ Real-Time PCR Thermal Cycler (Analytik Jena, Germany) with SYBR® Premix Ex TaqTM II (Tli RNaseH Plus) (Takara). The relative expression levels of mRNA were calculated with the ΔΔCT method^63^. The thermal cycling conditions were as follows: 95 °C for 3 min followed by 40 cycles of 95 °C for 15 s, 60 °C for 30 s and 72 °C for 30 s. Melting curve analysis from 65 to 95 °C was carried out after qPCR to ensure product specificity. The PCR amplification efficiency was analyzed by using different dilutions of the cDNA template. The standard curve for each gene was prepared according to the CT values at different cDNA concentrations. Primers with approximately 100% efficiency (R^2^ > 0.97) were used for qPCR (Table 2).

Statistical analyses

The CT values from the qPCRsoft for qTOWER^3^ (Analytik Jena) were analyzed in Microsoft Excel. IBM® SPSS® Statistics version 23 (SPSS, Inc., Armonk, NY, USA) was used for one-way analysis of variance (ANOVA), Duncan’s multiple range test (DMRT) and linear regression analysis. The residuals of CT values were evaluated by the difference between the actual value and the value calculated from linear regression for each gene (Fig. 1). The stabilities of the nine candidate reference genes were evaluated using geNorm version 3.5 (http://medgen.ugent.be/genorm/), NormFinder_0953 (http://www.mdl.dk), BestKeeper, RefFinder (https://blooge.cn/RefFinder/) and the comparative delta CT method. The raw CT values were transformed into 2^(-delta CT)^ for geNorm and NormFinder analysis. geNorm calculated the average expression stability value (M) for each gene, and compared the pairwise variation V to determine the optimal number of control genes for normalization. NormFinder software compared the estimated inter- and intra-group variances to identify the best gene with the lowest stability value. BestKeeper software used raw data, and ranked the stability by standard deviation (SD) and coefficient of variation (CV). RefFinder software integrated geNorm, NormFinder, BestKeeper and the comparative delta CT method, to rank the candidate reference genes with the geomean of the ranking values.

Supplementary Information

Supplementary Information.