All animal experimental procedures were approved by the Animal Care and Use Committee of Nankai University. In this study, samples of A. viridicyanea individuals were collected from Jizhou district (40.22′N, 117.50′E) in Tianjin, China on 19th May 2022 and were reared in an environmentally controlled growth chamber maintained at 25°C temperature with 60% relative humidity and a 16 h: 8 h (L:D) photoperiod, feeding on Geranium nepalense Sweet (Geraniaceae). Moreover, newly born eggs and hatched larvae of A. viridicyanea were obtained and kept in 90 mm diameter Petri dishes padded with moist filter papers, and then, hatched second larval instars were fed with normal host leaves or with those by antibiotic treatment until sexual maturation (10 days after emerging adults), as described in our previous report (Wei et al., 2020). In brief, 5 mL of antibiotic solution (50 μg/mL tetracycline, 200 μg/mL rifampicin, 100 μg/mL streptomycin = 1:2:4) was used to completely submerged leaves for 5 min. Then, a total of 50 normal (C group) and 50 antibiotic-reared individuals (A group), starved after 2 days, were used to gently dissect the whole gut tissues by clamping the head of living beetles with sterile forceps and removing the head. Finally, ten guts from the C or A group were surface-sterilized in 70% ethanol for 1 min, randomly pooled to build the transcriptome library as one replicate, immediately frozen in liquid nitrogen and stored at −80°C until subsequent total RNA extraction. In addition, the remaining guts after library construction were used for qRT-PCR analysis.

In this study, transcriptome sequencing was performed at Beijing Genomics Institution (BGI; Wuhan, China) Technology Co., Ltd. with three biological replicates of each treatment. Total RNA of all experimental guts was isolated using TRIzol reagent (Invitrogen, United States) according to the manufacturer’s protocols. In addition, RNA purity was assessed using the OD260/OD280 absorbance ratio between 1.8 and 2.0. RNA quality was visualized by 1% agarose gel electrophoresis staining, and RNA quantity was measured with an Agilent 2,100 Bioanalyzer (Agilent Technologies, United States). Finally, 5 μg of total RNA was removed from ribosomal RNAs (rRNAs) and used to prepare the transcriptome library construction.

Briefly, first-strand cDNA was produced with random hexamer primers and M-MuLV Reverse Transcriptase, and then the second strand was synthesized by adding RNaseH, dNTPs, DNA polymerase I, and buffer. Subsequently, the product was purified with AMPure XP beads, the cohesive ends of DNA were repaired using T4 DNA polymerase and Klenow DNA polymerase, and the product was further size-selected and ligated to an adaptor. Furthermore, we purified and amplified each product by AMPure XP beads and PCR (polymerase chain reaction), respectively. At the last step for preparation of the final sequencing library, magnetic beads were used to capture the amplified cDNA fragments, and the quality of all libraries was assessed using a Agilent BioAnalyzer 2,100 system (Agilent Technologies, United States). Later on, the diluted 1 ng/μL cDNA with an average insert size of 250–300 bp (base pair) was loaded for 150 bp paired-end sequencing using the MGI2000 platform.

The raw reads with adaptor contamination, low-quality, and undetermined data were removed using SOAPnuke software (version 1.5.2) to produce clean reads, and the Q20, Q30, and GC contents, representing effective sequencing values, were calculated. Furthermore, the quality of clean reads was estimated by FastQC (version 0.10.1), and then we mapped clean data to the genome of A. viridicyanea (data unpublished) with Hisat 2 (version 2.0.4) (Kim et al., 2019). The mapped reads of each group were assembled by StringTie (version 1.0.4), and finally, the transcriptome results were outputted (Pertea et al., 2015).

First, we applied Cufflinks software (version 2.2.1) to cover the transcripts with a clear chain orientation, and subsequently, the protein coding potentials were found by screening transcripts with no fewer than 2 exons and more than 200 bp in length. Moreover, CNCI software (CNCI_threshold >0 is mRNA, CNCI_threshold < 0 is lncRNA) was used to classify the coding transcripts and noncoding transcripts based on the spectrum of adjacent trinucleotides (Sun et al., 2013). Next, we used CPC (version 0.9-r2) (CPC_threshold >0 is mRNA, CPC_threshold < 0 is lncRNA) and txCdsPredict (txCdsPredict_threshold >500 is mRNA, txCdsPredict_threshold < 500 is lncRNA) software to assess the transcripts with coding potential according to the biological sequence characteristics (Kong et al., 2007). In addition, Pfam_Scan software was used to select the known protein family domains (Finn et al., 2010). In our study, the selected principle is that only transcripts identified by three of the above four methods were considered coding mRNAs or lncRNAs (Zhang et al., 2022b).

The expression levels of total mRNAs and lncRNAs were calculated and normalized to transcripts per million (TPM). In addition, RPKM (reads per kilobase of exon model per million mapped reads) was conducted to estimate their expression abundance in different samples, and the infrequently expressed transcripts with FPKM < 1 in all samples were filtered out. Furthermore, the correlation and principal component analysis (PCA) analyses based on RPKM using StringTie (version 1.0.4) software, and then the significantly differentially expressed mRNA and lncRNA were considered with absolute Log₂ | Fold Change | >1 and p-value < 0.001 by edgeR (version 3.14.0) software (Love et al., 2014).

First, to obtain annotation information on the functions of obtained mRNAs, all of them were searched against the Nr (NCBI nonredundant protein sequences), Nt (NCBI nucleotide), SwissProt (a manually annotated and reviewed protein sequence database), KOG (eukaryotic Ortholog Groups), and KEGG (Kyoto Encyclopedia of Genes and Genomes) databases using BLASTX (basic local alignment search tool) and BLASTN with an E-value < 10^–5. In addition, Blast2GO was used for GO (Gene Ontology) annotation with an E-value < 10^–5 from the Nr annotation results (Fu et al., 2021).

Moreover, the DELs and their target genes were functionally annotated by GO and KEGG databases. The GOseq R package was used to execute the GO enrichment, and KOBAS was used to explore the KEGG enrichment (Yu et al., 2012). The statistical significance between the A and C groups was considered the threshold to select significantly enriched GO terms and KEGG pathways (p-value < 0.01).

To explore the regulatory function of lncRNAs, cis/trans-regulatory algorithms were applied to predict the target genes. First, the regions 0–10 kb upstream or 0–20 kb downstream of cis-acting lncRNAs were screened to identify colocated target genes, and then the lncRNAs were considered trans-acting on the target genes once the binding energy was < −30 and exceeded the above region. Meanwhile, Pearson ≥0.6 and Spearman ≥0.6 correlation coefficients between the lncRNAs and target genes were used to predict the coexpressed target genes. In addition, the overlaps between the lncRNAs and target genes were classified into four categories: overlap, anti-overlap, incomplete, and anti-incomplete (Zhang et al., 2022b).

To identify candidate cDNA sequences of the Toll/Imd signaling pathway from sequencing data in A. viridicyanea, the published above genes from insects were used to search against the available information using TBLASTN and BLASTP databases with E-values < 0.001 as the threshold to find the corresponding gene members (Tong et al., 2015). Furthermore, the ORFs (open reading frames) of the putative cDNA sequences were translated using the ORF finder server, and the conserved domains of all the protein sequences were predicted by the SMART (simple modular architecture research tool) and TMHMM (version 2.0) servers. The signal peptide (SP) was examined using the program SignalP (version 5.0). In addition, for the protein sequence of each gene, basic information, such as molecular weight, length of protein, and isoelectronic point (PI), was identified with the ExPASy ProtParam server (Ren et al., 2021). Later on, the subcellular localization of the six Toll genes was predicted with WoLF PSORT (Zou et al., 2023). Lastly, MAFFT and MEGA (version X; Molecular Evolutionary Genetics Analysis) programs were applied to execute multiple sequence alignment of the amino acid sequences of the Toll/Imd signaling pathway genes from D. melanogaster, Anopheles gambiae, Bombyx mori, Tribolium castaneum, Diabrotica virgifera virgifera and A. viridicyanea, and to construct phylogenetic trees with the neighbor-joining (NJ) method (1,000 replicate bootstraps), respectively (Supplementary Table S2). Meanwhile, the results of phylogenetic analysis were optimized and shown using FigTree software (Tong et al., 2015). In this study, the obtained complete cDNA sequences of Toll/Imd signaling pathway genes were submitted to the GenBank database.

To our knowledge, lncRNAs, known as ceRNAs, can play vital roles in regulating specific target genes by sponging and sequestering miRNAs (Pasquinelli, 2012). Therefore, in order to better construct the ceRNA network based on the lncRNA-miRNA-mRNA axis, all the mRNAs of Toll/Imd signaling pathway were selected from the Section 2.7, and the principle of the selected related lncRNA and miRNA for building the ceRNA network is that the same expression trends of lncRNAs and mRNAs are opposite to the expression of miRNAs, and the p-values of lncRNA-miRNA (Pearson correlation coefficients ≥0.6) and lncRNA-miRNA (Pearson correlation coefficients ≥0.6) pairs were less than 0.05. Lastly, all these prediction was imported into Cytoscape software (version 3.9.1) to export the relevant ceRNA regulatory network (Otasek et al., 2019). Besides, the miRNA sequencing and prediction of the selected miRNA-mRNA pairs of Toll/Imd signaling pathway were in agreement with our previous report (data unpublished).

To integrate the transcriptomic and previous 16 S rRNA amplicon sequencing data (data unpublished), Pearson correlation coefficients were calculated for the integrated analysis of the above data. Concretely, the DEMs and DELs of Toll/Imd signaling pathway and the most abundant changes in gut bacterial microbiota from previous 16S rRNA amplicon sequencing were selected to produce correlation heatmaps using R packages (version 2.15.3). The level of statistical significance in this analysis was set at *p < 0.05, **p < 0.01 and ***p < 0.001.

To compare the expression levels between transcriptome sequencing and experimental analysis, fifteen DEMs and ten DELs of the Toll/Imd signaling pathway were used to perform quantitative real-time reverse transcription PCR (qRT-PCR). EF-1α was chosen as an internal control, and other primers were designed with Primer Premier 5 (Supplementary Table S3). Furthermore, the CFX96 real-time PCR detection system (Bio-Rad Laboratories, Hercules, CA) was used to carry out the 20 μL reaction system with PerfectStart Green qPCR SuperMix (Transgen, China) following the manufacturer’s instructions, including 10 μL of PerfectStart™ Green qRT-PCR SuperMix, 0.4 μL of each primer, 2 μL of cDNA template synthesized from total RNA in Section 2.1, and 7.2 μL of nucleotide-free water. The thermal cycling profile consisted of an initial denaturation at 95°C (30 s), followed by 40 cycles of 95°C for 15 s, 60°C for 30 s, and 72°C for 30 s. Finally, relative expression levels were calculated with the 2^−ΔΔCT method (Livak and Schmittgen, 2001), and the graphics of expression patterns were drawn by GraphPad Prism 9.0 (GraphPad Software Inc., United States).

In the present research, a total of six transcriptome libraries from the A. viridicyanea gut, namely, three control (C) and three antibiotic (A) treatment samples, were sequenced using the MGI2000 platform. Sequencing produced a total of 704,740,664 clean reads from 748,521,486 raw reads after removing low-quality, adaptor and N-containing reads. The mean Q30 and Q20 of each sample were 93.47% and 97.33%, respectively, and the mean GC content of the sequencing data was 47.82% (Supplementary Table S1). Furthermore, all these reads were used for transcript assembly and mapping with StringTie and Hisat 2, respectively, suggesting that different numbers of exons, RNA lengths (mRNA and lncRNA) and transcripts were obtained (Supplementary Figure S1). In addition, the clean reads were submitted to the NCBI SRA (sequence read archive) database under accession number PRJNA946198. All of the above sequencing parameters provided accurate evidence for further analyses.

Coding potential methods, including CPC, CNCI, Pfam, and txCdsPredict, were applied to screen the transcripts of mRNAs or lncRNAs. A total of 17,193 lncRNAs and 26,361 mRNAs (17,730 known mRNAs and 8,631 novel mRNAs) were obtained based on the prediction (Figure 1A, B). The lengths of lncRNAs and mRNAs, and the number of exons and transcripts are shown in Supplementary Figure S2.

At first, the FPKM values were calculated to estimate the expression patterns of all genes. In the C1 to C3 groups, the FPKM values were mainly less than 1, while the FPKM values were mainly between 1 and 10 in the A2 and A3 groups, except for the A1 group (Figure 1C). In our study, a heatmap was used to visualize all differentially expressed RNAs (Figure 1D). In addition, the correlation relationships (R ² values) represented the robustness of biological replicates and the reliability of RNA-seq data among the six groups (Figure 1E).

As shown in Figure 2, a total of 8,539 DELs were identified, including 8,120 upregulated and 419 downregulated DELs between the two experimental groups. In addition, 12,358 and 905 DEMs were upregulated and downregulated, respectively, from a total of 13,263 DEMs between the A and C groups. Of these DEMs, 8,432 and 3,926 were known and novel upregulated DEMs, and 650 and 255 were known and novel downregulated DEMs, respectively (Figure 2).

For target prediction of lncRNAs, lncRNA-mRNA interactions were screened based on the results of coexpression and colocalization. The results indicated that a total of 6,409 lncRNA-mRNA interaction pairs were cis/trans-acting from 6,035 lncRNAs targeting 6,421 mRNAs. In addition, a total of 4,086 lncRNA-mRNA pairs overlapped from 3,097 lncRNAs targeting 1985 mRNAs (Supplementary Figure S3). In addition, all DELs target a total of 6,425 DEMs in this research.

First, to obtain the basic annotation information of assembled mRNAs, all of them were compared against several databases. In detail, 25,594 mRNAs (97.09%) were annotated in at least one database, and 25,153 (95.42%), 13,072 (49.59%), 18,188 (69%), 21,278 (80.72%), 16,660 (63.20%) and 13,441 (50.99%) mRNAs were annotated against the Nr, Nt, SwissProt, KEGG, KOG, Pfam, and GO databases, respectively. Based on the Nr annotation, this transcriptome was mostly comparable to D. virgifera virgifera (50.65%), followed by Leptinotarsa decemlineata (14.4%), Anoplophora glabripennis (12.99%) and Callosobruchus maculatus (3.63%) (Figure 3A). Next, the GO database was used to comprehensively and widely investigate the basic functions of mRNAs, suggesting that a total of 13,441 mRNAs were classified into three GO categories, including biological process, cellular component, and molecular function (Figure 3B), such as 6,097 mRNAs in cellular process, 8,180 mRNAs in cellular anatomical entity, and 7,110 mRNAs in binding. Furthermore, the KEGG database was also used to systematically analyze gene functions. As shown in Figure 3A, C, total number of 18,188 mRNAs were assigned to six main categories of KEGG metabolic pathways: cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems, especially 1,638 mRNAs in the immune system.

Second, to ascertain the functional details of the DELs and target genes, GO and KEGG enrichment analyses were carried out. Based on the above GO categories, most DELs were distributed in cellular process, cellular anatomical entity, binding, catalytic activity, and so forth (Figure 4A). Likewise, 5,186 target genes were found in cellular anatomical entity, followed by 4,203 target genes in binding, 3,944 target genes in cellular processes, and 3,211 target genes in catalytic activity (Figure 4B). Furthermore, the top 20 significantly enriched pathways of DELs and target genes using KEGG enrichment indicated some immune signaling pathways, including the MAPK signaling pathway, TNF signaling pathway, B-cell receptor signaling pathway, IL-17 signaling pathway, and MAPK signaling pathway-fly (Figure 4C, D). In addition, we also summarized the various immune signaling pathways in the A. viridicyanea gut with different numbers of DELs and DEMs (Table 1). These collected data suggested that lncRNAs may contribute to the regulation of mRNA expression levels related to immune signaling pathways, which could help us to illustrate the immune responses of beetles under pathogen stimulation in the future.

In order to comprehensively understand one essential immune pathway in insects, we further focused on investigation of Toll/Imd pathway genes and their ceRNA network in A. viridicyanea based on sequencing data. Our observations suggested that six full-length cDNA sequences of Toll family genes were identified based on annotation information and online tools. The cDNA and protein characterization of all Tolls are shown in Table 2, and all sequences have been deposited in the NCBI database under accession numbers from OQ819300 to OQ819306. Concretely, all ORF lengths ranged from 3,972bp to 2,571 bp and encoded polypeptides 1,323 to 856 amino acids (aa) with molecular weights (MWs) of 149.45–98.26 kDa and theoretical isoelectric points (PIs) of 5.66–8.54, respectively. All Tolls contain different numbers of extracellular serial leucine-rich repeat (LRR) motifs, a transmembrane domain, and an intracellular cytoplasmic Toll/IL1 receptor (TIR) domain (Figure 5). Moreover, based on the classification of invertebrate Tolls described by Gerdol et al. (2017), six A. viridicyanea Tolls were divided into three subfamilies: four P-type TLRs (Toll-1 to Toll-4), one sPP-type TLR (Toll-5), and one sP-type TLR (Toll-6), whose subcellular location is the plasma membrane (Table 2).

Furthermore, to explore the signal members of the Toll/Imd pathway in A. viridicyanea, we also selected a total of nineteen extracellular and intracellular genes according to transcriptome sequencing. As shown in Table 3, six Spaetzle, three TAB2, two Relish, and one Myd88, Pelle, Tube, Cactus, Dorsal, Imd, Dredd and TAK1 genes were screened with complete coding DNA sequence (CDS) regions. Furthermore, we also characterized the above protein sequences in A. viridicyanea, including the length of the protein, molecular weight, isoelectronic point and functional domain, indicating that they possess typical domains for signal transduction or binding, such as the DEATH, Spaetzle, ANK, and RHD_DNA_bind domains. All these data showed that information on all Toll/Imd pathway genes with typical functional domains could expand our knowledge for comparative immunology between A. viridicyanea and other insects.

To evaluate the evolutionary relationships among the above family genes, the protein sequences from five insects were aligned to construct neighbor-joining trees (Supplementary Table S2). In line with our classification, the phylogeny result based on a phylogenetic tree with NJ method implied that the A. viridicyanea Tolls were well clustered together with corresponding subfamilies in insects, including P-type Tolls, sPP-type Tolls and sP-type Tolls (Figure 6A). All the pathway genes were also well clustered into each gene family (Figure 6B), which supports their orthologous relationships. Next, the immune recognition and signal transduction of all these genes will be validated in the A. viridicyanea gut.

Finally, to better understand the expression changes of Toll/Imd signaling pathway genes under antibiotic treatment, we calculated and summarized the profiles of DEMs and DELs of Toll/Imd pathway according to RNA-seq data. The results suggested that we found differential expression patterns of five Tolls, two Spaetzles and some pathway genes between the C and A groups. Most of these genes were significantly upregulated, while TAB2-1 expression significantly decreased during antibiotic rearing. Notably, the expression levels of all NF-κB-related genes were not significantly changed in this study (Figure 7A). In addition, the expression trends by qRT‒PCR were consistent with the differential expression of RNA sequencing, which verified the reliability of our sequencing data.

In this research, to investigate the miRNA sponge capacity of lncRNAs, eight Toll/Imd pathway genes, including Pelle, Dredd, Spaetzle-1, Spaetzle-3, Spaetzle-4, Spaetzle-5, TAK1 and TAB2-1, and their related miRNAs and lncRNAs were predicted and selected to construct a ceRNA network (Figure 8). In total, 76 miRNA‒mRNA pairs and 54 lncRNA‒miRNA pairs were obtained by target prediction. Of these lncRNAs and mRNAs, 10 DELs were involved in the regulation of 6 DEMs in this predicted lncRNA-mediated ceRNA network. In detail, Dredd, TAK1 and TAB2-1 are regulated by 7 miRNAs, and five miRNAs contribute to target Pelle. Meanwhile, 12 miRNAs were found to combine Spaetzle family genes. In addition, the novel miRNA-64-3p targets the majority of lncRNAs in the lncRNA-related ceRNA network of Toll/Imd signaling pathway. As shown in Figure 7B, a heatmap was used to display all above DELs with increasing expression between control and antibiotic groups, which brings together transcriptome sequencing and qRT-PCR validation with similar expression patterns.

We hypothesized that DELs and DEMs of the Toll/Imd pathway might be involved in response to alterations in the gut bacterial microbiota. Therefore, we detected their association relationships with Pearson correlation coefficients under antibiotic treatment. The results shown in Figure 9 suggest that, after antibiotic treatment, all selected changes in the gut bacterial microbiota were significantly positively or negatively correlated with different numbers of DEMs in the Toll/Imd signaling pathway. In detail, Toll-2 and Toll-3 changes were significantly positively affected by alterations in Bacteroides, Barnesiella, Prevotella, Parasutterella and Ruminococcus, whereas they significantly negatively influenced Toll-5 expression. In contrast, changes in Burkholderia, Eubacterium and Spiroplasma were only significantly correlated with Toll-4 and Toll-6 expression changes. Furthermore, changes in some gut bacterial microbiota, such as Bacteroides, Barnesiella, Prevotella, Parasutterella and Ruminococcus, were mainly significantly negatively involved in influencing some downstream factors, such as Myd88, Tube, Pelle, Imd, Dredd and Cactus, but Akkermansia alteration played the opposite role in affecting the expression of Pelle and Cactus. In addition, alterations in Eubacterium and Spiroplasma significantly positively or negatively affected the majority of the expression of DELs related to the Toll/Imd signaling pathway. However, Bacteroides, Prevotella, Ruminococcus and Akkermansia changes were only significantly positively or negatively associated with one DEL (Figure 9).