Integrated analysis of miRNA profiles and gut bacterial changes in Altica viridicyanea following antibiotic treatment

Abstract The gut bacteria involves in insect homeostasis by playing essential roles in host physiology, metabolism, innate immunity, and so forth. microRNAs (miRNAs) are endogenous small noncoding RNAs that posttranscriptionally regulate gene expression to affect immune or metabolic processes in insects. For several non‐model insects, the available knowledge on the relationship between changes in the gut bacteria and miRNA profiles is limited. In this study, we investigated the gut bacterial diversity, composition, and function from Altica viridicyanea feeding on normal‐ and antibiotic‐treated host plants using 16S rRNA amplicon sequencing; antibiotics have been shown to affect the body weight and development time in A. viridicyanea, suggesting that the gut bacteria of the normal sample were more diverse and abundant than those of the antibiotic‐fed group, and most of them were involved in various physical functions by enrichment analysis. Furthermore, we executed small RNA transcriptome sequencing using the two experimental groups to obtain numerous sRNAs, such as piRNAs, siRNAs, and known and novel miRNAs, by data mapping and quality control, and furthermore, a total of 224 miRNAs were identified as significantly differentially expressed miRNAs, of which some DEMs and their target genes participated in immune‐ and metabolism‐related pathways based on GO and KEGG annotation. Besides, regarding the regulatory roles of miRNA and target genes, a interaction network of DEM‐target gene pairs from eight immune‐ or metabolism‐related signaling pathways were constructed. Finally, we discovered that DEMs from above pathways were significantly positively or negatively correlated with gut bacterial alterations following antibiotic treatment. Collectively, the observations of this study expand our understanding of how the disturbance of gut bacteria affects miRNA profiles in A. viridicyanea and provide new valuable resources from extreme ranges for future studies on the adaptive evolution in insects.

microRNAs (miRNAs), serving as important regulators of gene expression, are a type of small noncoding RNA (sncRNA) with a length of approximately 18-24 nucleotides (nt) that are involved in a variety of biological processes, such as growth, development, apoptosis, and innate immunity, by promoting mRNA degradation or inhibiting translation to regulate the expression of target genes at the posttranscriptional level (Bartel, 2004;Zhao et al., 2022).
With the development of next-generation sequencing technology, it has become an essential tool to explore the miRNA profiles and their functional prediction through RNA sequencing technology, which provides an accessible opportunity to overall screen the differential expression of the genome at the transcriptional level (Wolf, 2013).To date, a growing number of studies have shown the functional roles of miRNAs in insects.For instance, appear to negatively regulate the Toll and Imd signaling pathway in Drosophila (Choi & Hyun, 2012;Lee & Hyun, 2014;Li et al., 2019;Li, Li, et al., 2017;Li, Xu, et al., 2017;Yao et al., 2023).In Bombyx mori (mulberry silkworm) larvae, miR-278-3p, interacting with IBP2 (insulin-related peptide binding protein 2), positively initiates the mRNA expression of cytoplasmic polyhedrosis virus (Wu, Qin, et al., 2016).It was shown that the JNK (Jun N-terminal kinase) pathway plays a key role in the immune system by regulating miR-184 in bacteria-infected pea aphids (Acyrthosiphon pisum; Ma et al., 2020).
Furthermore, the miR-n58/n174 cluster expression level was highly expressed in female white-backed planthopper, Sogatella furcifera, and meanwhile, the target gene, E3 ubiquitin-protein ligase RNF123, decreased in the females, as a sex-biased functional miRNA (Chang et al., 2018).In addition, it has been confirmed that Drosophila miR-92a/92b plays vital roles in the modulation of PDF neuronal excitability and maintenance of a neuroblast pool (Chen & Rosbash, 2017;Yuva-Aydemir et al., 2015).
The gut microbiota, such as bacteria, archaea, and eukaryotes, can produce beneficial components related to strengthening the immune system, defending against pathogens, and regulating host metabolism (Sekirov et al., 2010;Thursby & Juge, 2017).In herbivorous insects, symbiotic bacterial communities have likely evolved to favor species by detoxifying or resisting the toxins liberated in host guts, which is known as the gut microbial facilitation hypothesis (Hammer & Bowers, 2015).There is a report on the investigations of gut microbiome between the two Apriona species with host niche competition, suggesting that they were annotated into degrading plant toxic secondary compounds based on KEGG enrichment (Zhang, Wang, et al., 2022).Additionally, it is shown that the gut bacteria facilitate adaptation to crop rotation in the western corn rootworm (Chu et al., 2013), and they are involved in degrading defense chemicals consumed in various pests (Wei et al., 2020).Regarding the improving resistance against pathogens, in a hemipteran pest, Riptortus pedestris, the knockdown of the thanatin gene significantly increased the population of Burkholderia, a gut bacteria in the M4 crypt, in response to Escherichia coli K12 injection (Lee et al., 2022).The mortality of samples from a non-axenic Lymantria dispar asiatica larval gut was significantly higher than that of axenic larvae, and gut bacteria could be translocated from the gut to the hemocoel in the non-axenic group under Beauveria bassiana infection (Bai et al., 2023).To the best of knowledge, it is well known that the disturbance of microbiome homeostasis, known as "dysbiosis," resulted in an influence on physiological parameters within the immune systems or metabolic levels of host, owing to changes in the microbial community functionality (Egan & Gardiner, 2016;Legrand et al., 2020).
In recent year, the integrated analysis of multiomics data is available on the host-gut microbiota interactions in non-model dysbiosis insects.For example, antibiotic exposure disturbs the gut bacteria in Apis mellifera, and then influences host phenotypes, olfactory learning, and memory abilities using field or laboratory-generated bees combined with metagenomic sequencing, RNA sequencing, and proteome and metabolomics analysis (Zhang, Mu, et al., 2022).A recent study indicated that the gut bacteria of Plagiodera versicolora is an essential determinant of nutritional metabolism and immunity using comparative transcriptomic analysis and 16S rRNA amplicon sequencing of gut and body tissues in axenic and non-axenic larvae (Ma et al., 2021).Besides, Xie and colleagues observed miRNA and mRNA changes and constructed the miRNA-gene regulatory network under gut microbiota depletion using the abdomens of female Bactrocera dorsalis through RNA transcriptome and 16S rRNA amplicon sequencing (Xie et al., 2023).On the basis of the above observations, elucidating the disturbance of gut bacteria how to affect the target genes through regulation of miRNA expression will help us to deeply understand the adaptive evolution and maintaining homeostasis of insect populations.
It is easy to keep under laboratory conditions for investigating its physiological and behavioral changes, and its genomic data have been published (Xue et al., 2009(Xue et al., , 2016)).The gut microbiota is an essential determinant of adaptive evolution in A. viridicyanea, but there is little evidence on the interplay between gut bacteria and immune or metabolism-related miRNAs in A. viridicyanea combined with 16S rRNA amplicon and small RNA sequencing (sRNA-seq).In this study, we first used 16S rRNA amplicon sequencing by the Illumina HiSeq 2000 platform to detect the differences in the diversity and composition of gut bacterial communities in A. viridicyanea feeding on normal-or antibiotic-treated host plants.Furthermore, sRNA-seq analysis among above two experimental groups was performed on the A. viridicyanea guts to explore the characterization, expression patterns, and functional and regulatory roles of miRNAs.Lastly, a set of critical immune and metabolic differentially expressed miRNAs (DEMs) were selected, and the interaction relationships between DEMs and gut bacterial changes were executed based on bioinformatic prediction and Pearson correlation analysis.Taken together, | 3 of 16   REN et al.   this study can increase our understanding of the regulatory mechanisms of miRNAs and target genes in the insect gut and offer scientific data for future researches on the interaction between immunity or metabolism and gut bacteria in Coleoptera.

| Insect rearing and treatment
Approximately 40 adults of A. viridicyanea were collected from wild populations in Jizhou district (40.22′N, 117.50′E), Tianjin, P.R.China on May 19, 2022.All beetles were reared in a growth chamber held at 25°C with 60% relative humidity (RH) and long-day lighting conditions (16 h Light, 8 h Dark) and allowed to feed on its host plant G. nepalense, which was collected at the same location as insect samples.
Next, newly born eggs and hatched larvae of A. viridicyanea were collected and kept in dishes with moist filter papers under above rearing conditions.The molted 2nd-instar larvae were fed normal host leaves or host leaves completely submerged in antibiotic solution (50 μg/mL tetracycline, 200 μg/mL rifampicin, and 100 μg/mL streptomycin = 1:2:4) for 5 min, until 10 days after adult emerging, which were sexually mature adults (Wei et al., 2020).Furthermore, a total number of 80 normal-(control group, C) and 80 antibioticfed (antibiotic group, A) samples were starved for 2 days, and later, all whole gut tissues were gently dissected by clamping the head of live beetles with sterile forceps and removing the head.Finally, the gut tissues from 10 individuals of each experimental treatment were surface sterilized in 70% ethanol for 1 min and randomly pooled to construct the 16S rRNA or sRNA-seq libraries as one biological replicate.All samples were flash frozen in liquid nitrogen and stored at −80°C until subsequent experiments.In this research, all animal experimental procedures were performed and approved by the Animal Care and Use Committee of Nankai University.

| DNA extraction and 16S rRNA amplicon sequencing
A total of three biological replications from each treatment (n = 10 per replication) were used to extract DNA with a Genomic DNA extraction Kit (QIAamp) following the manufacturer's instructions.Furthermore, the amount, purity, and integrity of each DNA sample were estimated by a NanoDrop 2100 (Thermo) and electrophoresis of a 1.0% agarose gel, respectively.The results showed a single bright band on the gel of each DNA sample of approximately 20,000 bp (base pairs), and the absorbance ratio of A 260 /A 280 nm of all DNA samples ranging from 1.8 to 2.0 implied accessible DNA sample quality.Later on, the prepared DNA samples were diluted to 30 ng/μL with sterile water for 16S rRNA amplicon sequencing.
In brief, all samples were amplified with the same cycle numbers to ensure the accuracy and dependability of library construction using the barcoded primers V338F (5′-ACTCC TAC GGG AGG CAGCAG-3′) and V806R (5′-CGACT ACH VGG GTW TCTAAT-3′).Moreover, a QuantiFluorTM-ST blue fluorescence quantitative instrument was used to quantitatively detect the PCR (polymerase chain reaction) products.Next, the efficiency of the above PCR products was tested with 1.2% agarose gel electrophoresis, and then, they were purified with Agencourt AMPure XP beads (Beckman).The purified 16S rRNA amplicons were ligated with Illumina adapters (TruSeq DNA LT Sample Prep Kit), and the produced fragments with adapters were pooled in equimolar amounts.Finally, all the sequencing libraries from the above DNA products were built and sequenced with the Illumina HiSeq 2000 platform according to standard protocols at Beijing Genomics Institute (BGI) Technology Co., Ltd.

| Data analysis of 16S rRNA amplicon sequencing
In this study, we used FLASH (version 1.2.11) to filter the raw data for quality check and splice the paired-end sequences.Next, operational taxonomic units (OTUs) with sequence similarity were clustered through USEARCH (version 7.0.1090)software to discard chimeric sequences.RDP classifier (version 11.5) using the Bayesian algorithm conducted the taxonomic analysis based on the representative sequences of OTUs with a 97% similarity level using the I-sanger cloud data analysis platform.In addition, alpha diversity was estimated by QIIME (version 1.9.1) workflow (Caporaso et al., 2010), whose results presented the richness of gut bacterial communities by the Chao estimator.Besides, the ACE (abundance-based coverage estimator), Shannon and Simpson indices were calculated using Mothur (version 1.30.2) to reveal the diversity of gut bacterial community (Li, Wu, et al., 2021).The PCoA (principal coordinate analysis) based on Bray-Curtis dissimilarities was used to identify differences among gut bacterial communities.Collectively, the community heatmap was analyzed at the genus level, and Wilcox Test was applied to estimate the differences between the C and A groups.The functional annotation of bacterial OTUs was performed using Phylogenetic

Investigation of Communities by Reconstruction of Unobserved
States 2 (PICRUSt2), including COG (Clusters of orthologous groups of proteins), KEGG (Kyoto encyclopedia of genes and genomes), and Metacyc (Metabolic pathways from all domains of life) enrichment analyses (Douglas, 2015;Edgar, 2013).Finally, linear discriminant analysis effect size (LEfSe) was employed to identify biomarkers in different experimental groups (Segata et al., 2011).

| Total RNA extraction and sRNA-seq library construction
A total of three biological replications from each treatment (n = 10 per replication) were used to extract total RNA with TRIzol reagent (Invitrogen) according to the manufacturer's procedures.Moreover, the integrity, purity, and concentration of each RNA sample were estimated and assessed by 1.2% agarose gel electrophoresis and the OD 260 /OD 280 ratio by an Agilent 2100 Bioanalyzer system (Agilent Technologies), respectively.In addition, remaining RNA after library construction was used to qRT-PCR validation.
In the current study, for small RNA library construction, we built the sRNA-seq libraries from the qualified RNA using TruSeq Small RNA Sample Prep Kits according to the manufacturer's instructions.In brief, 1 μg of prepared total RNA was ligated to the 5′ and 3′ TruSeq adaptors and reversely transcribed by PCR.Next, the cDNA products were purified and selected from the 18 to 30 nt small RNA fragments by gel electrophoresis.The final constructed libraries were sequenced using 140-160 bp single-end reads on the DNBSEQ-G500 platform at Beijing Genomics Institution Technology Co., Ltd.

| Data assembly and filtration
Adaptor contamination, low-quality reads, and undetermined data were removed from the raw reads using SOAPnuke (version 1.5.0)software.Furthermore, the quality of clean reads was estimated by FastQC (version 0.10.1), and the information of Q20 and Q30, representing effective sequencing values, were statistically calculated.
Then, we used Bowtie (version 1.1.2) to map clean data to the reference genome of A. viridicyanea (Langmead et al., 2009).The mapped reads of each group were assembled by StringTie (version 1.3.0),and finally, the sRNA-seq data were produced (Pertea et al., 2015).
In addition, we used miRbase (version 22.0) and miRDeep2 (version 0.1.3)in a BLAST (basic local alignment search tool) search to identify known or novel miRNAs (Kozomara et al., 2019), and meanwhile, Piano and PHASIS (version 3.3) were used to predict piRNAs and siRNAs, respectively.

| Expression analysis of miRNAs
To characterize the expression patterns of identified miRNAs, the expression profiles were calculated and normalized to transcripts per million (TPM), and next, the correlation and principal component analysis (PCA) analyses were performed among all sRNA-seq libraries based on RPKM (Reads Per Kilobase of exon model per Million mapped reads) using StringTie (version 1.3.0)software and DESeq R (version 1.10.1)package, respectively (Pertea et al., 2015).The Bonferroni correction was used to adjust the p values in order to decrease FDR (false discovery rates, q values).The significantly differentially expressed miRNAs were those with absolute log 2 | fold change | >1 and p value <.001, calculated by edgeR package (version 3.14.0;Robinson et al., 2010).

| Target prediction and functional annotation
To further predict the target genes, RNAhybrid (version 0.1) and miRanda (version 3.3a) were applied to miRNA target prediction by identifying miRNA-binding sites with computational target prediction algorithms for both combined and overlapping genes (Enright et al., 2003;Lewis et al., 2005).Overlapping genes from the two algorithms were considered target genes.Moreover, to obtain annotation information on the functions of DEMs and target genes, we applied GO (Gene Ontology) and KEGG databases to functional enrichment analysis (Qu et al., 2023).Finally, we constructed a regulatory network of DEM-target gene axis from eight immune-and metabolism-related signaling pathways with Cytoscape software (version 3.9.1)based on the target and functional prediction (Otasek et al., 2019).

| Correlation analysis of DEMs and gut bacterial changes
In order to explore the relationships between sRNA-seq and 16S rRNA amplicon sequencing data, the Pearson correlation coefficients were calculated for the integrated analysis of the above sequencing results.Briefly, DEMs from eight signaling pathways and the top 10 changes in relative abundance of gut bacteria were selected to construct the heatmaps using R packages (version 2.15.3).Statistical significance was set at p values (*p < .05,**p < .01 and ***p < .001).
In addition, raw reads have been submitted to the NCBI database under accession number PRJNA946034 and PRJNA946198.

| Quantitative real-time PCR validation and statistical analysis
We randomly selected seven DEMs to validate the sRNA-seq results using quantitative real-time PCR (qRT-PCR).cDNA samples were synthesized from remanent total RNA from the Section 2.4 with TransScript® miRNA First-Strand cDNA Synthesis SuperMix (TransGen Biotech).Furthermore, qRT-PCR analysis was performed using TransScript® Green miRNA Two-Step RT-qPCR SuperMix (TransGen Biotech) on a CFX96 real-time PCR Detection System (Bio-Rad), and the following thermal cycling conditions were described as our previous study (Ren et al., 2023).The snRNA U6 was served as an internal reference gene, and other primers were designed by Primer Premier 5 (Table S9).All reactions were performed in triplicate, and the relative expression was calculated using the 2 −ΔΔc t method (Livak & Schmittgen, 2001).Finally, the column chart of relative expression of qRT-PCR and sRNA-Seq data was produced using GraphPad Prism (version 9.0; GraphPad Software Inc.).

| Gut bacteria in A. viridicyanea under antibiotic feeding
In the present study, we obtained a total of 386,148 clean reads with high quality from the normal-and antibiotic-fed A. viridicyanea guts after quality control using 16S V3-V4 rRNA amplicon sequencing.After alignment of these reads, a total of 488 OTUs were identified, and a Venn diagram of OTUs is shown in Figure 1a, suggesting that the two experimental groups shared 320 OTUs, and the C and A groups included 93 and 75 unique OTUs, respectively.Furthermore, the dominant taxa of the C group at the genus level was Salmonella, Raoultella, Barnesiella, and Bacteroides, while Barnesiella and Bacteroides were the dominant taxa in the A group, and the abundances of Salmonella and Raoultella were decreased (Figure 1c).Likewise, the gut bacteria of A. viridicyanea at the phylum level indicated that Proteobacteria and Bacteroidetes mainly existed in the C group, and the top 2 in the A group were Bacteroidetes and Firmicutes (Figure 1b).In addition, the most abundant genera and the hierarchical clusters of the gut bacteria are shown in the heatmap (Figure 1d).
We further calculated the alpha diversity of the gut bacteria between the C and A groups in A. viridicyanea.As shown in Figure 2, the Goods coverage index was more than 0.99 in the C and A groups, and the Chao, ACE, and Simpson indices in the C group were higher than those in the A group.Overall, the abundance and alpha diversity of gut bacteria between these two samples were altered, but the difference was not significant (p > .05).Moreover, the gut bacterial communities obtained from the six sequencing groups were significantly clustered into two different groups in the PCoA plot (Figure S1).Specifically, the taxa Bacteroideres, Bacteroidia, and Barnesiella were significantly enriched in the A group, and Delftia, Pedobacter, and Sphingobium were significantly enriched in the C group (Figure 2f).
Finally, to comprehensively illustrate the function of gut bacterial community in A. viridicyanea under normal and antibiotic treatments, we used the COG, KEGG, and Metacyc databases to classify and predict functions (Tables S1-S3).The data implied that translation, ribosomal structure, and biogenesis was the maximum functional classifications after antibiotic-fed treatment using COG enrichment analysis (Figure 3a).Carbohydrate metabolism, metabolism of cofactors and vitamins, and amino acid metabolism were the top 3 predicted pathways by KEGG annotation (Figure 3b), and nucleoside and nucleotide biosynthesis, amino acid biosynthesis, cofactor, prosthetic group, electron carrier, vitamin biosynthesis, and fatty acid and lipid biosynthesis were the main metabolic functions predicted by Metacyc annotation (Figure 3c).All these predicted functions may represent the most critical functions of the gut bacteria, and they played an essential role in host homeostasis.group) and 94,961,924 (A group) raw reads, which were filtered to 98,807,924 and 92,356,655 clean reads, respectively, with average Q20 and Q30 scores of more than 98.97% and 96.87% (Table 1).The resulting lengths of clean reads are shown in Figure S2, ranging from 18 to 32 nt, suggesting that the miRNA length range is 18-24 nt, piRNA length range is 26-32 nt, and siRNA length range is 20-25 nt.

| Data analysis of small RNA sequencing
In this research, we further focused on the identification, expression, and function of miRNAs.

| Identification of known and novel miRNAs
The results revealed that the mapping ratio of clean reads exceeded 86.49% in each sRNA-seq group by mapping onto the A. viridicyanea genome (Table S4).Moreover, all clean reads were compared to known sRNAs with various databases, indicating that a total of 122,766 and 168,803 reads were mapped to sRNAs in the A and C groups, respectively, showing that 39 sRNAs were regarded as rRNAs and tRNAs identified.Herein, we identified a total of 2220 miRNAs, including 44 known miRNAs and 2176 novel miRNAs, with U (uridine) bias at the first nucleotide position and G (guanine) at the least frequent nucleotide (data not shown).

| Expression patterns of miRNAs
Principal component analysis (PCA) revealed strong aggregation of all sRNA-seq samples (data not shown), and besides, the R 2 values (Pearson correlation coefficients) were 0.891-1.000among the six sRNA-seq groups based on FPKM values, representing the robustness of the biological replicates and reliability of the sRNA-seq data (Figure 4a).Furthermore, as shown in Figure 4b, the TPM values of  most miRNAs were less than 10 among the six sRNA-seq groups.The gut tissues of antibiotic-fed individuals were used as experimental groups, and normally reared samples were used as controls, suggesting that a total of 224 miRNAs were screened as significantly differentially expressed miRNAs under the condition of p value <.001 and log 2 | fold change | >1, of which 127 miRNAs were upregulated and 97 miRNAs were downregulated (Figure 4c and Figure S3).
Finally, seven DEMs were randomly selected to qRT-PCR validation, revealing that the expression levels of DEMs by qRT-PCR validation were consistent with the sRNA-seq data (Figure S4).

| Target prediction and functional annotation
To explore the functional roles of miRNAs and target genes, we iden- targeting 10,594 genes.The results demonstrated that target genes were involved in a wide range of biological functions, including biological process (BP), cellular component (CC), and molecular function (MF) by GO annotation.In detail, 47 GO terms were assigned to target genes; of these terms, 22 terms belonged to biological process, 14 terms belonged to cellular component, and 11 terms belonged to molecular function, including binding, followed by cell, cell part, cellular process, catalytic activity, and so forth (Figure 5a and Table S5).
As shown in Figure 5b, DEMs were mainly enriched in binding, cellular process, cell, cell part, and catalytic activity.Notably, 12 DEMs were enriched in immune system process, 136 DEMs were enriched in response to stimulus, and 81 DEMs were enriched in multicellular organismal process.
To clarify the functional pathways of DEMs and target genes, the KEGG enrichment results suggested that a total of 5152 target genes were annotated to 44 KEGG pathway subcategories from six categories, including cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems (Figure 6, Table S6).The dominant pathways of target genes were global and overview maps, signal transduction, cancers: overview, and infectious diseases: viral.2).In summary, these results indicated that the annotation information of DEMs and target genes could lay the foundation for further in-depth studies to characterize their roles in insect homeostasis.

| Interaction network of selected DEM and their target genes
Based on the functional enrichments and target prediction, we ran- diamond) was found to regulate only 15 target genes in the JAK/ STAT signaling pathway but not those in metabolism-related pathways (Figure 7b).This finding suggested that some DEMs could play multifunctional roles in the regulation of physical activities and provided new basis for further exploring the regulatory mechanisms based on the competing endogenous RNA (ceRNA) networks in A. viridicyanea.

| Correlation analysis between selected DEMs and gut microbial changes
To determine the correlation relationships between two sequencing data, we selected the top 10 abundant changes in gut bacteria and 10 DEMs from each signaling pathway (the selected pathways are similar to Section 3.6) using Pearson correlation analysis.At first, detailed information of randomly selected DEMs from eight pathways for correlation analysis is presented in Table S7, containing various upregulated or downregulated DEMs in each pathway.Moreover, the data is shown in Table S8 and Figure

| DISCUSS ION
There is accumulating evidence that the homeostasis between the insects and gut microbiota community was maintained within diverse biological processes (Kim & Lee, 2014).The gut microbiota contributes to effect on the insect's physiological traits at both the local (gut) and systemic (body) levels (Ma et al., 2021).Without this knowledge, the molecular mechanisms of host-microbiota interactions in non-model insects could not be deeply elucidated.
Here, we speculated that the dysbiosis of gut bacterial microbiota might affect miRNA expression to further regulate target genes using normal-and antibiotic-fed A. viridicyanea, and it has been revealed that the gut bacterial alteration in beetles by antibiotic treatment resulted in a significant decrease in the body weight of newly emerged adults and prolonged development time (Wei et al., 2020).
Collectively, our data mainly investigated the relationships between changes in the gut bacteria and DEMs and laid the foundation for the potential mechanisms for gut microbiota-insect interactions in Coleoptera.
Diverse gut bacteria in insects can significantly impact the host homeostasis (Douglas, 2015).For example, the larvae of A. swainsoni possess highly organized gut bacteria in four gut segments that benefit its host, including adapting to the relatively harsh niche conditions and playing a major role in host growth (Zhang, Zhuang, et al., 2022), so we aimed to get more detailed knowledge about the variation in diversity, composition and function of gut bacteria from normal-and antibiotic-reared A. viridicyanea.Based on 16S rRNA amplicon sequencing, the number of sequences per sample was far greater than 60,000, and numerous OTUs were identified.The alpha diversity of gut bacteria from the normal samples was more abundant than that in the antibiotic-fed groups, such as the Chao, ACE, and Simpson indices, but there were no significant differences among the two experimental groups, which is in agreement with our previous report (Wei et al., 2020).The composition of the gut bacterial community of A. viridicyanea was further investigated at the phylum level, with Firmicutes, Bacteroidetes, and Proteobacteria being the most dominant in the two experimental groups, which can affect a variety of aspects of host biology, such as nutritional absorption, insecticide resistance, and the degradation of defense chemicals (Ley et al., 2006;Salem et al., 2017;Zheng et al., 2023).Furthermore, the organization of gut bacterial community at the genus level was also observed, implying that Salmonella and Raoultella, as well as Bacteroides and Barnesiella, were most prevalent in the C and A samples, respectively.Mounting data indicated that they may be involved in inflammation, aging, and carbohydrate and shortchain fatty acid metabolism ( | 13 of 16 REN et al.Li, Yang, et al., 2021;Zhao et al., 2023).Finally, we annotated the functional information of gut bacterial community from normal-and antibiotic-fed A. viridicyanea when the results were compared with three databases, mainly including organic metabolism and biosynthesis.Likewise, it is well known that KEGG significantly enriched in the gut bacteria of A. germari and A. swainsoni contained biosynthesis and metabolism-related pathways (Zhang, Wang, et al., 2022).All NAs with approximately 18-24 nt and a dominant length of 22 nt, which was similar to previous reports (Chang et al., 2018;Wu, Jiang, et al., 2016;Wu, Qin, et al., 2016;Yang et al., 2022), piRNA with 26-32 nt in length, and siRNA with 20-25 nt in length.Later on, we identified a total of 2220 miRNAs, containing 44 known miRNAs and 2176 novel miRNAs, among which 224 DEMs were calculated between the A and C groups, consisting of 127 upregulated and 97 downregulated DEMs.Moreover, target gene prediction was performed with RNAhybrid and miRanda, illustrating that ten thousand target genes were obtained in this study.
The GO database provided useful information for understanding gene functions and specific processes, indicating 923 target genes in binding, followed by 698 target genes in cell and cell part, 682 target genes in cellular process, and 653 target genes in catalytic activity.
Furthermore, total number of 203 DEMs were classified in binding, followed by 199 DEMs in cellular process, 196 DEMs in cell and cell part, and 191 DEMs in catalytic activity and metabolic process, which is consistent with the results in other insects (Wu et al., 2013;Wu, Jiang, et al., 2016;Wu, Qin, et al., 2016).Moreover, we also found that target genes and DEMs involved in various signaling pathways, including immune processes and metabolism, according to KEGG enrichment.As we known, insects can activate the humoral response, for example, the production of antimicrobial peptides by Toll/Imd (immune deficiency) and downstream signaling pathways, lysozymes, and rapidly activated phenoloxidase (PO) cascade-mediated melanization constitute the humoral immune response of insects, and besides, they also initiate the cellular response to phagocytose bacteria and encapsulate parasites (Lemaitre & Hoffmann, 2007;Myllymaki et al., 2014;Sheehan et al., 2018).In Drosophila, the immune responses against gram-positive bacteria and fungi mainly depend on the Toll pathway, whereas the Imd pathway has more function in defense against gram-negative bacteria by the NF-κB signaling pathway (Choe et al., 2005;Lemaitre et al., 1995;Leulier et al., 2002).In addition, there are two parallel pathways in the Drosophila gut, the DUOX pathway, and Imd pathway, which play indispensable roles in controlling the immune response and maintaining homeostasis in fruit flies (Bosco-Drayon et al., 2012;Ha et al., 2005).
In previous data, in vivo experiments demonstrated a requirement for Ras/MAPK signaling in restricting innate immune responses in hemocytes, fat bodies, and adult gut stem cells in Drosophila (Ragab et al., 2011).The PI3K/Akt signaling pathway was significantly enriched with KEGG enrichment analysis, and upregulated genes from this pathway were identified in Holotrichia parallela, and moreover, silencing PI3K resulted in increasing susceptibility to entomopathogenic nematode and Bacillus thuringiensis stimulation and accelerating symbiotic bacterial multiplication (Li, Wu, et al., 2023).Here, our study provided more valuable data for the regulation between miRNA and immune pathways in non-model insects.In addition to immune function, miRNAs have impact on the host metabolic systems, for instance, glutathione metabolism and metabolism of xenobiotics by cytochrome P450 contain cytochrome P450 family genes and glutamate and glutamine genes, which play important roles in the transcriptional regulation of hypoxia, maintaining normal immune system function and antioxidant effects in animals (Saetan et al., 2020).On the basis of the above observations, to better understand the regulatory relationships of DEMs and target genes, a co-ex-

F
I G U R E 1 (a) The Venn diagram shows the number of shared and specific OTUs from the C and A groups.Average relative abundances of gut bacteria at the phylum (b) and genus (c) levels (Taxa with abundances <1% are included in "others").(d) Heatmap clustering analysis of the gut bacteria at the genus level.Next, linear discriminant analysis (LDA) effect size (LEfSe) analysis was used to detect discriminating taxa between the C and A groups (p < .05 and LDA > 2), indicating that the gut bacteria of A. viridicyanea in the antibiotic-fed group was represented by some biomarkers.
In this study, six sRNA-seq libraries from the normal-and antibiotic-fed A. viridicyanea guts generated a total of 102,216,241 (C F I G U R E 2 The alpha diversity of gut bacteria in A. viridicyanea from the C and A groups is estimated based on the number of OTUs, with coverage (a), Chao (b), ACE (c), Shannon (d), and Simpson (e) indices using Wilcox Test.(f) The taxa with significantly different abundances of gut bacteria of A. viridicyanea based on linear discriminate analysis (LDA) score from the C and A groups.

F
The functional classifications of gut bacteria in A. viridicyanea using COG (a), KEGG (b) and MetaCyc (c) databases.All the x-axes represent the two experimental groups, and the y-axes represent the relative abundance.Various colors indicate the term names in different databases.TA B L E 1 Quality assessment of sRNA-seq data.
2014 miRNAs targeting 11,395 genes, including 224 DEMs F I G U R E 4 (a) Correlation heatmap of sRNA-seq among the six sRNA-seq groups.(b) The expression levels of all miRNAs among the six sRNA-seq groups.(c) Volcano plot of DEMs between the C and A groups.Each point in this figure represents one miRNA.The x-axis indicates the gene expression with fold change (antibiotics vs. control), and the y-axis indicates the p values adjusted by the false discovery rate (FDR; q value) levels.F I G U R E 5 (a) GO enrichment of target genes.Each bar represents the gene number in each GO term in the category of BP (biological process), CC (cellular component) and MF (molecular function).The x-axis indicates the number of genes.(b) GO enrichment of DEMs.A larger circle indicates a greater number of enriched DEMs.The x-axis indicates the comparison group.All the y-axes are the GO functional classifications.

Furthermore
, we identified and summarized numerous DEMs from various immune and metabolism-related pathways, such as the Toll/ Imd signaling pathway, Toll-like receptor signaling pathway, NODlike receptor signaling pathway, RIG-I-like receptor signaling pathway, NF-κB signaling pathway, fatty acid metabolism, metabolism of xenobiotics by cytochrome P450, and so forth (Table
8, demonstrating that Bacteroides, Barnesiella, Prevotella, Parasutterella, and Ruminococcus are not significantly positively correlated with DEMs in the Toll/Imd and NF-κB signaling pathways, and they are less significantly associated with DEMs than those in the JAK/STAT and MAPK-fly signaling pathways.Similarly, Bacteroides, Barnesiella, Prevotella, Parasutterella, and Ruminococcus were also not significantly positively associated with DEMs in four metabolism-related pathways.Overall, the number of DEMs in eight pathways significantly negatively correlated with changes in the gut bacteria was larger than that with significantly positive correlation.Additionally, some DEMs, such as novel-miRNA-1066-5p, novel-miRNA-879-5p, novel-miRNA-981-5p, novel-miRNA-872-5p, and novel-miRNA-970-5p, were conversely significantly association with different gut bacterial changes in same pathway.
these results suggested that changes in the gut bacterial community by antibiotic treatment may affect the physiological biochemical action of host insects.MiRNAs are one kind of sncRNAs that are approximately 22 nt in length, and they can influence immune responses, development, cell proliferation, and apoptosis in animals by base pairing with the target mRNA 3′-untranslated regions (3′-UTR) to trigger mRNA degradation or translational suppression(Gebert & MacRae, 2019;Li & Chen, 2021).In the current research, we constructed the six sR-NA-seq libraries from normal-and antibiotic-fed A. viridicyanea guts to explore the miRNA profiles through next-generation sequencing technology.SRNA-seq data produced massive clean reads with high quality, and numerous sRNAs were screened, including miR- pression network was constructed according to enrichment data and target prediction, showing that numerous genes could be potentially regulated by different DEMs in eight signaling pathways, of note, the novel miRNA-196-5p might play a certain function in innate immune regulation in the A. viridicyanea gut, and further, their function should be deeply investigated in the future.To our knowledge, the gut bacteria of insects can play important roles in homeostasis to improve host fitness by supplying essential amino acids and nutrients and affecting host immunity against invading pathogens, for example, some hemipteran insects have to abandon some innate immune systems, including the AMPs and Imd pathways, F I G U R E 8 Correlation analyses between the gut bacteria and DEMs of immune or metabolism-related pathways from the C and A groups.Each column represents different gut bacteria, and each row represents various DEMs.Red and blue colors represent negative or positive correlations, respectively.Ellipses sizes are proportional to correlation value, and correlation significance is shown by *** (p < .001),** (p < .01)and * (p < .05).(a) Toll/Imd signaling pathway, (b) JAK/STAT signaling pathway, (c) MAPK-fly signaling pathway, (d) NF-κB signaling pathway, (e) metabolism of xenobiotics by cytochrome P450, (f) carbon metabolism, (g) fatty acid metabolism, and (h) glutathione metabolism.

Sample ID Raw reads Total clean reads Low quality reads rate (%) Invalid adapter reads rate (%) Read length too small rate (%) Clean reads Q20 (%)
Summary of DEMs in immune and metabolism-related pathways.