A susceptible strain (S) of SSB was sourced from the Institute of Plant Protection, Chinese Academy of Agricultural Sciences (Haidian District, Beijing), provided by Dr. Lanzhi Han and reared on an artificial diet (Shuijin et al., 2017) under constant laboratory conditions (temperature: 27°C ± 1°C, relative humidity: 70%–80%, and photoperiod: 16:8) without insecticides for over 50 generations. The S strain was screened under chlorantraniliprole treatment for over 30 generations, and the resistance ratio (RR) of the R strain was 110.4 compared to the S strain, as previously reported (Sun et al., 2018). The SSB of both S and R strains were reared on an artificial diet with or without chlorantraniliprole treatment.

Total RNA was extracted from three to four 4th instar SSB larvae using the TRIzol reagent kit (Life Technologies, United States), following the manufacturer’s instructions. RNA quality was first checked by 1% (w/v) agarose gel electrophoresis, followed by a NanoDrop spectrophotometer (RNA model, California, United States). Lastly, Qubit 2.0 accurately quantified the RNA concentration. We finished the library construction and RNA-seq (Novogene, China). The total RNA of six samples (three independent biological replicates for S and R strains, respectively) was purified and qualified first, followed by the removal of rRNA using the Ribo-Zero kit (Epicentre, United States), the random disturbance of RNA in fragmentation buffer (NEB, United States), and the cDNA synthesis with 6-bp random hexamers. The purified double-stranded cDNA was repaired and adapter-ligated. Then, the cDNA library was enriched by PCR after being purified. Sequencing was performed on the NovaSeq 6000 platform.

To isolate SSB lncRNAs from the sequencing dataset, a computational pipeline was constructed, which was referred to report with minor modifications (Liu et al., 2017). At the very beginning, Trimmomatic software was used to filter low-quality reads (Bolger et al., 2014). The raw reads from the six libraries were then mapped to the SSB genome via TopHat (Trapnell et al., 2009). In detail, the reads from each library were first mapped to the scaffolds, and the junction outputs from every RNA-seq dataset were combined together to generate a “pooled junction set,” which was further used to map all of the reads from different RNA-seq datasets to the scaffolds via TopHat. This step supplied a junction set for Cufflinks (Trapnell et al., 2012). Then, the six datasets were combined into a whole transcriptome by Cuffcompare, according to the annotated genomic information. The transcripts that are longer than 200 nT and have more than two exons were reserved, and 24,267 transcripts were obtained in this step. At this point, the candidate protein-coding genes were excluded after blasting to the NR database (e < 001). Next, transcripts >300 nT ORF were removed by getorf software (http://emboss.sourceforge.net/apps/cvs/emboss/apps/getorf.html). Then, the protein-coding transcripts were predicted by Coding Potential Assessment Tool (CPAT), and other transcripts were used as the template to search in the Pfam database with HMMER software (Finn et al., 2015). The transcripts that have no conserved domain/motif potential were reserved, while other non-coding RNAs except lncRNAs, the known tRNAs, snoRNAs, snRNAs, and rRNAs were all removed by the Infernal and BLASTN search (Zhao et al., 2016), generating the final lncRNA sets.

Aligning lncRNAs with the SSB genome (Ma et al., 2020), we analyzed the lncRNA exon–intron structure by Geneious (Kearse et al., 2012) and also checked its distribution among the scaffolds.

The transcript enrichment of the identified lncRNAs was examined by counting reads and normalized by Cuffdiff software, using the t-test to display the significance of DE. The lncRNA expression level was detected with FPKM (fragments per kilobase of transcript per million fragments mapped) as an indicator. Q-value is the FDR-adjusted p-value. A lncRNA which meets these criteria will be defined as specifically expressed: 1) the value of |log2 fold-change|≥1; 2) p-value <.01; and 3) q-value<.01. The fold change refers to the ratio of expression amount between two strains. The lncRNAs obtained by the differential expression analysis are displayed in the form of a volcano plot.

Sample preparation, RNA extraction, and quality control are as described previously. Three reactions during the cDNA synthesis are forward (F) primer with reverse transcriptase (RT), reverse (R) primer with RT, and F + R primers without RT. The primers for RT-PCR (Supplementary Table S1) were designed using http://www.idtdna.com/Scitools/and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China).

The extraction and purification of total RNA were described previously. cDNA was obtained by using the PrimeScript^®RT Master Mix kit (TaKaRa, Japan) with 500 ng of total RNA as the template. qRT-PCR was carried out with a SYBR Premix Ex Taq kit (TaKaRa, Japan) on an ABI 7300 real-time PCR system (Thermo Fisher Scientific, United States) as described as follows: 30 s denaturation at 95°C, 40 cycles of 5 s at 95°C, and 31 s at 60°C, 15 s at 95°C, 15 s at 60°C, and 15 s at 95°C again. Actin A1 was used as an internal reference. Data were analyzed by the 2^−ΔΔCt method (Pfaffl 2001). The primers used are shown in Supplementary Table S2.

Samples of different development stages were collected including the egg, 1st, 2nd, 3rd, 4th, 5th, and 6th stages, pupae, and male and female moths. For each stage, 3–4 individuals were collected as one biological replicate, and four biological replicates were performed. These samples were then used to conduct RNA extraction, cDNA collection, and qRT-PCR for the lincRNA MSTRG.7482.1.

The chemically synthesized small interference RNA (siRNA) came from GenePharma Co., Ltd. (Shanghai, China) (see Supplementary Table S2 to find its sequence).The HPLC-purified double-stranded siRNAs were dissolved in diethylpyrocarbonate-treated water (Milli-Q-grade) to make a 4 mg/ml solution. Then, 1 μl (4 μg) of siRNA was injected into fourth-instar larvae with a micro-needle, pulling the needles from glass capillaries which have 1.0-mm outer diameter and 50-mm inner diameter by using a micropipette puller (Model P-87, Sutter Instruments Co., CA), keeping needles still for 30 s at the injection point to avoid siRNA leakage. The shuffled siRNAs were taken as negative control (Supplementary Table S2). About 20 to 30 specimens were examined for each treatment, with triplicate experiments. The fourth larvae were collected 24 and 48 h post-micro-injection to determine the efficiency of RNAi.

For RNAi assay, the larvae were removed onto an artificial diet without any pesticide treatment for 24 h post-injection. The larvae that died or whose movement was blocked due to mechanical injury were discarded. The remaining specimens were removed onto an artificial diet with chlorantraniliprole treatment, according to the LC₅₀ dose in S or R strains for another 96 h, respectively. Lastly, these specimens were removed onto an artificial diet without treatment. The mortality was calculated day by day from the fourth day post-treatment with chlorantraniliprole.

For the location of lncRNA and protein-coding genes on the genome of SSB, information on the chromosome length, lncRNA/gene density, and location of the lncRNA/gene to be located was collected. The R package RIdeogram was used to display the aforementioned information. Regarding the qRT-PCR result, Tukey’s test and one-way ANOVA statistical analysis were performed to compare the different gene expression levels between S and R strains.

Strand-specific RNA-seq was finished via NovaSeq 6000, Illumina. Raw data on RNA-seq were submitted to the public database: National Genomics Data Center (https://ngdc.cncb.ac.cn/). The assigned accession of the submission is CRA009124. Then, 43,490,487 to 55,131,982 clean reads were obtained from a total of 94.73 G clean data (Supplementary Table S3). Mapping these reads to an updated version genome of SSB (Ma et al., 2020), the highest alignment rate of 91.85% was yielded. An average of 121,998 transcripts was assembled in each sample from aligned reads. Finally, a total of 40,573 merged transcripts were determined (Supplementary Table S3). Twenty-two transcripts shorter than 200 nT were first removed from 40,573 merged transcripts (Figure 1A). Then, 16,284 transcripts which contained single exons were also discarded. Up to 24,267 transcripts with multiple exons were filtered by CPAT, Swiss-Prot, and Pfam databases to exclude the possibility of protein-coding genes. Finally, a total of 3470 lncRNAs were discovered and characterized in SSB after removing small non-coding RNAs via NONCODE and Rfam databases. These lncRNAs were further divided into four categories, namely, intronic lncRNAs, lincRNAs, sense lncRNAs, and antisense lncRNAs, according to their relative positions and orientations on the genome (Figure 1B). Most lncRNAs (1879, 54.12%) were located in intergenic regions, the second-ranked sense lncRNAs (853, 24.57%) shared an overlap with the exon region of protein-coding genes, and 14.23% and 7.09% lncRNAs were antisense lncRNAs (493) and intronic lncRNAs (245) (Figure 1B). Among 3470 SSB lncRNAs, 2454 of which contained two exons, accounting for 70.68% (Figure 2A), and only 2.13% of lncRNAs contained seven or more exons. In contrast, only 13.27% protein-coding genes harbored two exons. In addition, the protein-coding genes with three exons accounted for the largest proportion (17.16%). The exons whose length was shorter than 200 bp made up for more than 70% in protein-coding genes, but there was less than 40% in exons shorter than 200 bp in lncRNAs (Figure 2B). As shown in Figure 2C, the average length of the SSB lncRNA transcript and the protein-coding genes were 1049 and 199 bp, respectively.

Up to 3356 lncRNAs were gathered in S or R strains after removing 116 lncRNAs with FPKM=0 in all six libraries. There were 1343 lncRNAs dysregulated with |log2 fold-change|>1. Also, there were 798 and 545 lncRNAs highly expressed in R and S strains, respectively (Figure 3A). Moreover, a total of 632 differentially expressed (DE) lncRNAs met strict screening criteria, that is, p-value < 01, p-adj < 01, and |log2 fold-change|>1 (Figure 3B).There were 297 and 335 lncRNAs upregulated in S and R strains. Among these lncRNAs, 365, 31, 141, and 95 belonged to lincRNAs, intronic lncRNAs, sense lncRNAs, and antisense lncRNAs, respectively.

In the detoxification enzyme gene family, P450s, carboxyl/cholinesterases/esterases, and ABC genes were related to the chlorantraniliprole resistance in SSB (Wang et al., 2015; Xu et al., 2019; Peng et al., 2021). Also, a 1000-kb range between the center of an lncRNA gene and a neighboring protein-coding gene transcription start site (TSS) was defined as a regulatory zone by the Genomic Regions Enrichment of Annotations Tool (GREAT) (McLean et al., 2010). To study the potential function of lncRNAs in chlorantraniliprole resistance, both the detoxification enzyme genes and DE lncRNAs were located in the genome of SSB. Sixty-six P450 genes were located on 20 different chromosomes (Figure 4A; Supplementary Table S4). Among the 38 lncRNAs adjacent to these P450 genes, the lincRNA, sense lncRNA, antisense lncRNA, and intronic lncRNA were 24, 7, 5, and 2, respectively. On chromosome10, the lincRNA MSTRG.2776.2 was upregulated in the S strain as 8.45 log2 fold change (log2FC) as in the R strain, and it was 944, 528 bp apart from CYP9A12. For esterases, there were 67 carboxyl/choline esterases located on 16 chromosomes (Figure 4B; Supplementary Table S5), and a total of 27 lncRNAs were located nearby. On chromosome 9, the lincRNA MSTRG.25727.1 was upregulated in the R strain as 6.54 log2FC as in the S strain, and CsuEst35 antennal esterase was 887, 454 bp apart from it. Forty-three of the 47 ABC transporter genes previously reported were successfully located in the genome of SSB. They are distributed on 16 different chromosomes (Figure 4C; Supplementary Table S6). At the same time, 51 lncRNAs were adjacent to them, including 31 lincRNAs, 12 sense lncRNAs, 6 antisense lncRNAs, and 2 intronic lncRNAs. ABCB4, D3, and G12 were located on chromosome 8, and the expression of D3 was upregulated and that of G12 was downregulated in the R strain (Peng et al., 2021). A total of 18 lncRNAs were dysregulated in the R strain (Figure 4C; Supplementary Table S6) as lincRNAs MSTRG.25315.3 and MSTRG.25316.8 were upregulated.

Considering lincRNAs, also known as long intergenic non-coding RNAs, accounted for 54.10% of the number of lncRNAs (Figure 1B, 1879 out of 3470), our focused analyses revealed there were 167 and 198 lincRNAs differentially highly expressed in S and R strains, respectively. Those lincRNAs whose value of |log2FC| ranked in the top 20 were chosen in the following analysis (Tables 1, 2). The log2FC of lincRNA MSTRG.6875.40 was 12.44, ranking it the highest in the top in the S strain (Table 1), whereas MSTRG.25315.3 was 15.70, ranked first in the R strain (Table 2). The log2FC values of differentially expressed MSTRG.22483.6 and MSTRG.13055.1 were 3.82 and 6.16, as they were listed last. These 40 lincRNAs were distributed on 17 different chromosomes (Chr1, 2, 4, 6, 8, 9, 10, 12, 13, 15, 16, 17, 18, 21, 26, 27, and 28) or scaffolds (scaffold 789 and scaffold 83). The length of them varied widely, ranging from 231 to 11,712 bp. The exon numbers varied from 2 to 7. For example, there were seven, six, and five exons in MSTRG.25315.3, MSTRG.11012.32 (Table 2), and MSTRG.6875.40 (Table 1), respectively.

Protein-coding genes within 1000 kb adjacent to lincRNA were collected. MSTRG.25315.3 was the top differentially upregulated gene in the R strain (Table 2, 15.70 of log2FC) with 83 coding genes located. Heat shock protein 70–2 (GenBank accession: AGR84224.1) was only 1,770 bp close to MSTRG.25315.3 (Figure 5). The lincRNA MSTRG.25727.1 was 337 bp long, but it harbored four exons and was elevated to a log2FC value of 6.54 in the R strain compared to the S strain. Among 20 mRNA genes nearby, antennal esterase was 887, 454 bp adjacent to MSTRG.25727.1. In addition, sodium channel protein was 429, 858 bp adjacent to MSTRG.7200.3 (Table 2).

With regards to the top 20 highly expressed lincRNAs in the S strain, MSTRG.2776.2 was 961, 794 bp adjacent to cytochrome 9A20 (annotated from Bombyx mori, GenBank accession: BAI47532.1). Sixty-four protein-coding genes were close to MSTRG.6875.40 within 1000 kb, and UDP-glucuronosyl transferase 2B31-like was 990, 115 bp apart from it. MSTRG.20713.2 was upregulated in the S strain (Table 1, 7.84 of log2FC), and it was 494, 273 bp away from 4A1-like gene, which is a family member of solute carrier organic anion transporters, encoding membrane proteins that are involved in solute (charged and uncharged organic molecules, and inorganic ions) transportation. They are likely to be involved in pesticide metabolism and transport.

We randomly selected 20 lincRNA candidates shown in Tables 1, 2 to be validated by strand-specific RT-PCR. Results showed 13 lincRNAs were strand-specific, that is, five came out forward stranded (MSTRG.25723.1, 3932.19, 3932.24, 7482.1, and 6875.40) and eight were reverse stranded (MSTRG.8464.1, 25727.1, 11012.32, 29804.8, 29804.10, 22483.6, 25315.3, and 25316.8) (Figure 6). Of the 13 strand-specific lincRNAs, only four were in the same direction as the originally labeled strand (Tables 1, 2, MSTRG.25316.8, 29804.8, 25727.1, and 29804.10). The other seven lincRNAs can be amplified by both primers (Figure 6).

To verify the RNA-seq results, the relative expression levels of three randomly selected DE lincRNAs in the R strain (MSTRG.25315.3, MSTRG.25316.8, and MSTRG.17788.5) and five lincRNAs in the S strain (MSTRG.3932.19, MSTRG.3932.15, MSTRG.7482.1, MSTRG.29805.1, and MSTRG.22483.6) were examined by qRT-PCR (Figure 7). Consistent results were observed between the qRT-PCR results of most of the selected strains and the sequencing data except for MSTRG.17788.5 and MSTRG.7482.1. The two lincRNAs displayed the opposite trend. The expression of MSTRG.17788.5 was 9.23 times the log2FC value in the R strain as that in the S strain (Table 2), but it turned out 0.69 times in the R strain as that in the S strain. On the contrary, the expression of MSTRG.7482.1 was first 4.92 of log2FC in the S strain as that in the R strain (Table 1), but it turned out 35.86 times in the R strain confirmed by qRT-PCR.

MSTRG.7482.1 had the highest differential expression folder change among the three lincRNAs, which is highly expressed in the R strain (Figure 7). We first used qRT-PCR to examine the expression of MSTRG.7482.1 in different stages of development, like eggs, each instar larvae, pupae, and male and female moths (Supplementary Figure S1). Altogether, MSTRG.7482.1 expressed highly from eggs until the sixth-instar larvae. Then, it declined sharply in the pupal stage and was expressed to the lowest in the female moth. The expression profile suggested MSTRG.7482.1 functions mainly in larvae, which is, the feeding stage. The siRNA was designed from position 231 based on the sequence of MSTRG.7482.1, and we successfully knocked down MSTRG.7482.1 gene by injecting siRNA in the fourth-instar larvae (Figure 8). The expression of MSTRG.7482.1 was knocked to about 40% compared to control in both 24 and 48 h after injection.

To detect the effect of the downregulated MSTRG.7482.1 gene on chlorantraniliprole sensitivity, the fourth-instar larvae injected with siRNA-231 were removed onto an artificial diet with chlorantraniliprole treatment, according to the LC₅₀ dose in the S or R strain for 96 h, respectively. In the S strain, the RNAi group showed 70.09% average mortality, which is significantly higher than that of control (mortality: 38.62%) four days after pesticide treatment (Figure 9A). On day 8 post-chlorantraniliprole treatment, both mortality of siRNA-231 and the control group declined and it turned out closer (39.36% vs. 26.75%) (Figure 9B). In the R strain, the average mortality of the siRNA-231 group was higher than that of the control at 4 and 8 days after pesticide treatment; there was no statistically significant difference (Figures 9C, D).