Unravelling the gut bacteriome of Ips (Coleoptera: Curculionidae: Scolytinae): identifying core bacterial assemblage and their ecological relevance

Bark beetles often serve as forest damaging agents, causing landscape-level mortality. Understanding the biology and ecology of beetles are important for both, gathering knowledge about important forest insects and forest protection. Knowledge about the bark beetle gut-associated bacteria is one of the crucial yet surprisingly neglected areas of research with European tree-killing bark beetles. Hence, in this study, we survey the gut bacteriome from five Ips and one non-Ips bark beetles from Scolytinae. Results reveal 69 core bacterial genera among five Ips beetles that may perform conserved functions within the bark beetle holobiont. The most abundant bacterial genera from different bark beetle gut include Erwinia, Sodalis, Serratia, Tyzzerella, Raoultella, Rahnella, Wolbachia, Spiroplasma, Vibrio, and Pseudoxanthomonas. Notable differences in gut-associated bacterial community richness and diversity among the beetle species are observed. Furthermore, the impact of sampling location on the overall bark beetle gut bacterial community assemblage is also documented, which warrants further investigations. Nevertheless, our data expanded the current knowledge about core gut bacterial communities in Ips bark beetles and their putative function such as cellulose degradation, nitrogen fixation, detoxification of defensive plant compounds, and inhibition of pathogens, which could serve as a basis for further metatranscriptomics and metaproteomics investigations.

Apart from providing commodities such as food and wood, forests perform fundamental ecological services such as biodiversity preservation and climatic regulation 1 . Forest health in Europe is under severe threat due to global climatic change. The trade-off in tree carbon investment between primary and secondary metabolites is observed under abiotic stress conditions such as drought 2 . On top of impacting the tree defence physiology adversely, repeated years of drought also trigger more recurrent and devastating bark beetle outbreaks leading to large scale forest decline 3,4 . However, this is not true for healthy trees as conifers can produce entomotoxic defensive chemicals, such as monoterpenes, in addition to releasing anti-aggregation compound such as 4-allylanisole to defend themselves against bark beetle attacks [5][6][7][8][9][10] . Conifer defence primarily composed of terpenoid resins made up of mono-, sesqui, and diterpenes among which monoterpenes are already reported entomotoxic in higher concentrations to bark beetles 11,12 . Moreover, conifer tissue feeding is nutritionally limiting for bark beetles because of the low concentration of phosphorus, nitrogen, sterols, and vitamins 13 . Successful bark beetle colonization demands strategies to overcome the noxious effect of plant secondary metabolites such as avoidance of toxin ingestion by modification of feeding behaviour; plant defence manipulation by converting the highly toxic metabolites to low toxic ones; enhanced excretion of ingested toxins; sequestration of the toxin; target-site mutation or metabolic degradation of toxin compounds 14 . Using aggregation pheromone mediated mass attack, some bark beetles (European spruce bark beetle, Ips typographus) can surmount host defence as the conifers unable to accumulate sufficient defence response to beetles arrived en masse 15  Gut bacterial diversity in bark beetles. OTU abundance. The bacterial sequences obtained from the gut tissue of five Ips and one non-Ips bark beetles were clustered into a total of 3049 OTUs at 97% similarity cutoff limit (Supplementary excel 3). The rarefaction curve and estimation of good's-coverage indicating the completeness of sampling (> 99%) represent the whole bacterial diversity for each of the five Ips and one non-Ips bark beetle species that were covered (Supplementary Figure 1, Table 1). The sequences thus obtained were assigned to 43 bacterial phyla. Among these 43 bacterial phyla, the predominance of 5 bacterial phyla, namely Proteobacteria, Firmicutes, Actinobacteria, Bacteroidetes, and Tenericutes ( Fig. 1) was represented using GraPhlAn. The relative abundance of the top 100 genera identified in this study was represented in the evolutionary tree as well as the top 10 genera was further denoted in the taxonomic tree (Supplementary Figure 2A,B). In particular, considering the top 10 phyla, the relative abundance of Proteobacteria was higher in all the beetle species (IT-96%, ID-79.6%, SX-95.4%, IAC-54.1%, IC-55.4% and PP-96.2%). However, IAC gut documented a relatively high abundance of Firmicutes (31.6%) and Bacteroidetes (4.7%) whereas Tenericutes (18.4%) and Actinobacteria (19.4%) were prevalent in IC ( Fig. 2A  2) (P < 0.001). Moreover, notable differences in gut-associated bacterial community diversity among the beetle species were also documented in the present study. IAC feeding on pine trees (Shannon index-4.8 ± 0.6; Simpson index-0.8 ± 0.1) showed significant differences (P < 0.05) compared to the spruce feeding Ips beetles (IT, ID) indicating the influence of different host feeding on shaping up the gut-associated bacteriome ( Fig. 3 and Table1). However, the extent of bacterial diversity may vary from beetle species to species. Besides, the influence of host feeding, the significant difference (P < 0.01) between the pine feeding beetles IAC (Shannon index-4.8 ± 0.6; Simpson index-0.8 ± 0.1) and SX (Shannon index-1.8 ± 0.2; Simpson index-0.6) may also depict the environmental impact on the gut bacteriome. Furthermore, the number of observed bacterial species was higher in IAC (680.2 ± 38.4) compare to the other Ips bark beetles collected from R site (IT, ID, IC, and SX) ( Table 1) indicating the plausible location-specific influence on gut bacterial richness but this needs further experimental validation. Considering the Ips bark beetles feeding on spruce (IT, ID), 274 core OTUs were observed, whereas the pine feeding bark beetles from two different sites shared 213 OTUs in common (Fig. 4A,B and Supplementary excel 5,6). Furthermore, the present study documented the occurrence of 126 core OTUs shared among all five Ips species (IT, ID, IC, SX and IAC) (Fig. 4C) that were assigned to 44 families and 69 genera. The shared community www.nature.com/scientificreports/  . Apart from the core bacterial community present in the bark beetle gut, a consortium of different unique OTUs was also detected. However, it is important to note that unique OTUs do not always refer to the unique bacterial species; nonetheless, the diversity of the unique bacterial species follows a similar trend in the present study (Supplementary excels 5,6,7).
β-diversity. Beta (β) diversity evaluating the differences in microbial communities between the bark beetles feeding on the different coniferous hosts was represented by the box plot and heatmap based on Unweighted UniFrac and Weighted UniFrac distances between samples (Supplementary Figures 4, 5A). The hierarchical clustering based on Unweighted UniFrac distance matrix (UPGMA-Unweighted pair group method with arithmetic mean) illustrates the influence of environment on gut bacterial diversity where all Ips beetle species collected from R site (IT, ID, IC, and SX) were clustered together in one clade (Supplementary Figure 5B). Similarly, Non-Metric Multi-Dimensional Scaling (NMDS) result also showed the clustering of the beetles together from the same site regardless of feeding on different conifer hosts suggesting the likely impact of sampling location on bark beetle gut microbial community (Fig. 5). Conversely, Weighted UniFrac UPGMA tree clustered all three spruce feeding beetles together irrespective of their location of sampling (Supplementary Figure 5C). Hence, no generalization can be derived at this point regarding the degree of influence sharing between the host plant and the location of bark beetle sampling in shaping up the gut bacteriome. The extent of the top 12 bacterial species-specific variations at the genus level (Metastat analysis) revealed the predominance of Erwinia in Ips bark beetles collected from R-site (IT, ID, IC, and SX) while Sodalis is highly abundant in non-Ips PP collected from K-site. Considering the spruce feeding Ips beetles, the high abundance of Rahnella and Acinetobacter was detected in IT, whereas Listeria and Tyzzerella were abundant in ID. However, among the pine feeding bark beetles, Tyzzerella, Pseudoxanthomonas, Faecalibacterium, Subdoligranulum, Nocardioides, and Acinetobacter were most abundant in IAC whereas Serratia and Rahnella dominate in the gut of SX (Fig. 6). Moreover, t-test analysis performed to compare species-specific variation between the spruce feeding www.nature.com/scientificreports/ Ips bark beetles IT and ID from the same collection site revealed a significant difference in abundance of Serratia and Pantoea (Supplementary Figure 6). Furthermore, ANOSIM analysis evaluating the variation between the bacterial community among the six different bark beetles was significantly higher than the variation within the group (Supplementary Table 1). Our data revealed that the bacterial communities present in five Ips and one non-Ips bark beetle species are significantly different. However, MRPP and ADONIS analysis determining significant differences of the overall gut bacteriome among the bark beetles showed no such differences between IT and ID, both collected from R site and feeding on the spruce (Supplementary Tables 1, 2). In contrast, AMOVA analysis documented no significant differences between IT and SX (Supplementary Table 3) indicating marginally significant differences need to be interpreted carefully. Linear discriminant analysis effect size (LEfSe) indicated the differentially abundant bacterial species between the bark beetles. The histogram of the LDA (Linear discriminant analysis) scores represented the predominance of the gut bacterial community of beetles feeding on different conifer hosts. Bark beetle species were characterized by a preponderance of some of the following significantly abundant genera (LDA score [log10] > 4) such as Erwinia, Rahnella, Raoultella, Pseudomonas, Stenotrophomonas in IT; Serratia in ID; Actinobacteria, Spiroplasma, Enterococcus in IC; Serratia, Rahnella, Rickettsia in SX; Tyzzerella, Vibrio, Pseudoxanthomonas, Curtobacterium, Streptococcus, Faecalibacterium in IAC while Sodalis, Wolbachia, Pantoea, Enterobacter in PP (non-Ips). (Supplementary Figure 7A). Considering all six bark beetles together, the presence of 1 to 7 bacterial biomarkers per beetle species was observed (Supplementary Figure 7B). Furthermore, within the spruce Ips feeding beetles (IT, ID) members from Proteobacteria (class-α-Proteobacteria and γ-Proteobacteria) and Firmicutes (class-Clostridia and Bacilli) served as predominant biomarkers. Similarly, bacterial species from the phylum Bacteroidetes (class-Sphingobacteriia) was also prevalent in pine feeding beetle along with Firmicutes  www.nature.com/scientificreports/ functional diversity of the bacterial communities' present in the beetles where the contribution of each OTU is associated with a given gene function. The relative abundance of the following gene functions such as carbohydrate metabolism, cell growth and death, immune-related function, lipid metabolism, metabolism of Terpenoids and Polyketides, biodegradation and metabolism of xenobiotics, biosynthesis of other secondary metabolites, amino acid metabolism, transport and catabolism, signalling of molecules and their interaction, revealed at level-2 KEGGs Orthologs (KOs), were highly represented in IAC. The genes involved in metabolic diseases, replication, and repair, translation, nucleotide metabolism, immune response, energy metabolism was dominant in larch feeding beetle IC (Fig. 8A). However, considering the relative abundance of the top 10 gene functions, no significant difference was observed among different Ips bark beetles (Fig. 8B). Such profile of highly abundant gene function such as membrane transport, carbohydrate metabolism, amino acid metabolism, replication and repair, cellular processes and signalling, translation, metabolism of cofactors and vitamins and nucleotide metabolism might be attributed to core bacterial communities in Ips, which may contribute to the fundamental functions in bark beetle physiology. Comparing with Ips and non-Ips bark beetles, the genes involved in metabolic diseases, degradation, environmental adaptation, metabolism of cofactors and vitamins, transcription, glycan biosynthesis, and metabolism and immune function were relatively less abundant within Ips bark beetles.

Discussion
Bark beetles serve as a keystone species in forest ecosystems by nutrient cycling of mature and wind-felled trees. However, these beetles become a severe pest for production forests during the outbreak phase, causing considerable economic as well as environmental damage. The present study characterizes the gut-associated bacterial communities of five Ips and a non-Ips bark beetles from the Scolytinae subfamily feeding on different conifer hosts. The study reveals the core bacterial communities in the gut of two spruce feeding, two pine feeding, and one larch feeding Ips bark beetles collected from two forest locations within Czech Republic and their putative functional relevance in bark beetle holobiont. In our study, the bacterial community structure reflected by the α-diversity indices illustrates lower bacterial diversity and richness in Ips typographus (IT), Ips sexdentatus (SX), while Ips acuminatus (IAC) exhibited highest richness and diversity among all Ips bark beetles. ID and IC displayed similar bacterial diversity and richness that is somewhat in between IT, SX and IAC (Fig. 3). The higher diversity of bacteriome in IAC may also be co-related with their fungi feeding behaviour during larval stages 53 . High bacterial diversity may be obligatory to support such feeding habit. Furthermore, IAC being a less aggressive beetle compare to Ips sexdentatus considering their ophiostomatoid fungi arsenal 50 , may also need a wide array of bacterial species to cope with the pine defensive compounds. However, considering the attack behaviour in the forest IAC can attack healthier trees than SX 54 , www.nature.com/scientificreports/ which is considered as typical secondary pest 55 . Alternatively, Ips typographus is considered as aggressive beetle species that execute a mass attack to overwhelm the plant defence mechanism and thereby can sustain well with low bacterial diversity in their gut. However, such inferences would need further experimental validation where both aggressive and less aggressive bark beetle species need to be collected from a particular location, feeding on the same host tree to nullify the influence of microenvironments and the differences in feeding habitat. Besides, identification of the host tree microbiome would also be informative in this context. Bark beetle gut bacterial assemblage is a non-random process 56 . However, several factors such as beetle sampling location, population characteristics (epidemic or endemic), competition for survival and resource acquisition, host quality or availability, and trophic interactions within the microhabitat can directly or indirectly influence the bacterial community at a given time in the bark beetle gut. Despite all these sources of variabilities, our results document the presence of a core bacteriome representing nine phyla identified to 44 bacterial families and 69 bacterial genera across all five Ips bark beetles (Fig. 4C). It can be presumed that the obtained core microbiome involved in conserved function in bark beetles. The data represents the predominance of Proteobacteria followed by Firmicutes, Actinobacteria, and Bacteroidetes that are consistent with the previous studies done on other bark beetles from the same subfamily 57 . To colonize and survive under the tree bark, the bark beetles largely depend on their gut bacterial species to digest the complex celluloses, hemicelluloses, and lignin to provide nutrition to the beetles 31,33 . In the current study, bacterial genera, namely Erwinia, Serratia, Rahnella, Raoultella, and Pantoea belonging to Enterobacteriaceae family, frequently observed in the core bacteriome in Ips bark beetles may contribute to such conserved function to survive in the similar habitat under the bark 57 . Besides, Pseudomonas, Stenotrophomonas, Pseudoxanthomonas, Novosphingobium, Sphingomonas and species from Ruminococcaceae family also present among the core microbiome are reported to degrade cellulose and other plant cell wall polysaccharides such as starch, xylan, and lignin, and the breakdown products derived may www.nature.com/scientificreports/ provide nutrition to beetles 31,33,34,58,59 . Furthermore, Acinetobacter detected in our study demonstrate esterase and lipolytic activity 60 . Such abilities of the bacterial species to hydrolyze the phospholipids, triglycerides, and fatty acids present in the phloem and resin of the tree 61 may contribute to energy reserve in the fat bodies of bark beetles that can be utilized later during beetle development as well as host colonization and reproduction 62,63 .
Handling the nitrogen limitation in the diet is a communal challenge for all bark beetles feeding on a nitrogenlimiting phloem diet. Bark beetles often rely on their gut symbiotic bacteria capable of nitrogen acquisition for mitigating the crisis 64 . Nitrogen-fixing bacteria belonging to genera Rahnella, Stenotrophomonas, Pantoea, Raoultella, Pseudomonas, Burkholderia, and Bradyrhizobium, which are present as the core bacteriome may be abetting Ips bark beetles out of the crisis [65][66][67][68] . Moreover, Pseudomonas and Rahnella can recycle the uric acid from the beetle faeces into assimilable nitrogen, aiding bark beetle sustenance 65,66 . For successful tree colonization, the bark beetles also face challenges to combat the host defensive compounds that are released in high amounts soon after bark beetle attack 69 . Though such compounds at low levels enable the bark beetles in finding their suitable host 5 , the same host defence compounds in high concentration (i.e., terpenoids) are entomotoxic. To outcompete the plant defence mechanisms, some bark beetles follow different strategies such as overwhelming the www.nature.com/scientificreports/ tree defence by a mass attack or by surviving in the hostile environment with assistance from their symbionts 16 . Pseudomonas, Serratia, Rahnella, Erwinia, and Pantoea identified as core members of Ips bark beetle gut in the present study may help in the deprivation of toxic monoterpenes in addition to degradation of the plant cell wall polysaccharides and nitrogen fixation in those beetles 35,36,41 . Moreover, bacterial symbionts are often implicated in providing vitamins and amino acids to insect host 20,70,71 . Some bark beetles ingest fungi that supplement the beetles with vitamins, amino acids, and sterols 27,29 . Interestingly, bacterial species Pseudomonas is reported to produce a variety of antifungal antibiotics to prevent the growth of pathogenic fungi in Dendroctonus genus 72,73 , whereas Stenotrophomonas and Pantoea are documented to function against the pathogen in red turpentine beetles 66 . Additionally, Clostridiaceae observed in the core bacterial community is also reported to produce antimicrobial molecules, but their function is yet to be discovered 74 . Furthermore, Pseudomonas, Serratia, Rahnella, and Erwinia that are detected in the present study is reported to be associated with the production of bark beetle anti-aggregation pheromone, verbenone 40 . Plant defensive compound such as α-pinene serves as a precursor of verbenone, a common bark beetle anti-aggregation pheromone 75,76 . The biochemical pathway of verbenone production from cis-verbenol is recently documented to be mediated by enzymes originated from facultative anaerobes in D. valens gut (red turpentine beetle) under the frass-simulated and gut-simulated oxygen concentration environment 77 . Our present findings also endorse such possibilities; however, that requires further investigation. Our data propose that the Ips bark beetles share a persistent consortium of core bacterial communities that may contribute to the fundamental metabolic pathways such as cellulose degradation, nitrogen fixation, detoxification of defensive plant compounds, and inhibition of pathogens. However, the variation in the gut bacterial communities among the bark beetles, observed from the β-diversity, maybe due to demographic factors that comprise the discrepancies in feeding habitat and other environmental characteristics besides the prevailing difference at the species level. The environmental influence (i.e., location of sampling) on the bark beetle ecology plays a vital role in shaping the gut bacteriome. Interestingly, the Ips bark beetles collected from the R site (IC, SX, IT, and ID) are clustered together (Fig. 5) and show similar variations overall in their gut bacterial communities when compared to other Ips bark beetle from L-site (IAC) and non-Ips bark beetle (PP) collected from K-site. Such findings could be explained by the fact that several bacterial species are also acquired from feeding on the different coniferous hosts in the same forest environment. However, such conclusions need further experimental corroborations.
However, it is interesting to note that besides the core bacterial communities, the Ips bark beetles also host different bacterial genera in the gut, which could be assimilated during feeding on to the phloem tissue. The occurrence of Spiroplama belonging to Tenericutes is significantly higher in larch feeding beetle, IC (18.4%) compared to other beetles (< 0.5%) but completely absent in IAC. Previous studies reported Spiroplasma as an insect symbiont, but its role in the gut is yet to be identified 78 . Furthermore, Pseudomonas, Serratia, and Pseudoxanthomonas, present as core gut members in IC, are also detected in the cuticle and galleries of other larch feeding beetle, Dendroctonus simplex (eastern larch beetle) 79 . The presence of Actinobacteria in the gut of the IC may produce antimicrobial compounds. It is documented that the presence of Streptomyces belonging www.nature.com/scientificreports/ to Actinobacteria phylum is producing antifungal agents against the growth of Ophiostoma minus in Dendroctonus frontalis Zimmermann (southern pine beetle) 80 . Based on our knowledge, it is the first record for Sodalis to be present in the gut of bark beetles. However, Sodalis is only present in IAC (12.4%) and absent in other Ips beetles. Though the role of Sodalis in bark beetles is not yet understood completely, hitherto, it is described as a mutualistic endosymbiont in tsetse flies and heteropteran insects with extensive metabolic capabilities 81,82 . In addition to the above-mentioned bacterial species, the predominance of Wolbachia and Rickettsia are considered as a symbiotic pathogen associated with several beetle species for inducing cytoplasmic incompatibility, resulting in reproductive distortions and hence, may serve as candidates for future management of bark beetles 83,84 . Noteworthy, each beetle species hosts a set of specific bacterial genus as biomarkers, and some of the bacterial biomarkers congregate in the functional profile of the bark beetles evident from the PICRUSt analysis (Fig. 8). Some of the biomarkers detected in the present study have the potential to play a vital fundamental role in bark beetle physiology hence can be the targets for further investigations. For instance, the presence of bacterial species belonging to the families Lachnospiraceae, Cellulomonadaceae, Ruminococcaceae, Chitinophagaceae is associated with degradation of complex plant polysaccharides by their cellulolytic enzymes 58,85 . In addition to the carbohydrate metabolism ability, some members of Enterobacterials and Pseudomonadales are also capable of nitrogen fixation, terpenoid degradation and metabolism of xenobiotic compounds 35,36,41,66 . The production of antimicrobial compounds to inhibit the growth of pathogens and contributing to the immune response of the bark beetles are attributed to the presence of bacterial biomarkers such as Pseudomonadales 72,73 . Moreover, the top 10 physiological functions contributed by gut bacterial community seems to be conserved among all Ips beetles, and that is perhaps associated with the core bacterial communities (Fig. 8B). However, such functional aids from the gut bacteria to the host beetles need to be experimentally validated.
In conclusion, it is a pioneering study surveying the gut-associated (endomicrobiome) bacterial community in five Ips bark beetle species inhabited in Czech forests. Sixty-nine core bacterial genera are documented and explored for their putative ecological functions. Our data advocates plausible conserved functions for the core gut-bacteriome in different Ips beetles irrespective of ecological and demographic variabilities among them. Hence, they can serve as targets for future downstream functional studies (i.e., metatranscriptomics and metaproteomics). It is essential to mention here that the in-depth knowledge about the core gut bacterial communities and their fundamental function in Ips bark beetle adaptation to survive in a hostile environment could expand our knowledge about the important forest nutrient recycler and may unleash the unexpected potential for keeping bark beetle population in the non-epidemic stage.  88 . To assimilate six biological replicates per beetle species representing the beetle population in the area, more than 120 living and healthy beetles were collected and pooled from infested logs in each locality belonging to 8 or more different trees and subsequently shock frozen under liquid nitrogen for future use. It is noteworthy to mention here that due to the pooling of beetles during collection, the individual colony specific variability of the beetle gut microbiome, which was not the primary study goal, could not be estimated.

Materials methods
Randomly selected bark beetles were surface sterilized and dissected using dissecting scopes in the biosafety cabinet under sterile conditions 89 . Only guts that were free of apparent nematode infections were used for downstream DNA extraction. Subsequently, gut tissues (8 to 10 guts pooled per replicate) were homogenized, and microbial DNA was extracted using PureLink Microbiome DNA Purification Kit (Invitrogen, Thermo Fisher Scientific, US) following manufactures protocol. Purified DNA was quantified using Qubit 2.0 Fluorometer (Thermo Scientific) using Qubit dsDNA HS Assay kit (Invitrogen, Thermo Fisher Scientific, US), and the integrity was evaluated by 1% agarose gel electrophoresis. Lastly, extracted high-quality DNA from all beetle species (six biological replicates per species) was sent for high-throughput amplicon sequencing (Novogene Company, China). Amplicon sequencing. Amplicon sequencing was performed following the pre-optimized protocol. Precisely, DNA was diluted to 1 ng/μL using DNase free water, and bacterial 16S rRNA genes of distinct regions (V3-V4) were amplified using a specific primer (341F-806R) 90 with the barcode. All PCR reactions were performed with Phusion High-Fidelity PCR Master Mix (New England Biolabs). PCR products were visualized on 2% agarose gel, and samples with clear amplification between 450 and 480 bp were selected for further experiments. PCR products were purified with the Qiagen Gel Extraction Kit (Qiagen, Germany). Sequencing libraries were generated using NEBNext Ultra DNA Library Pre-Kit for Illumina, and index codes were added. The library quantity and quality were measured on the Qubit 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. Lastly, the library was sequenced on a Novaseq 6000 platform, and 250 bp paired-end reads were generated. www.nature.com/scientificreports/ Data analysis (α-diversity/ β-diversity and core microbiome/functional analysis/statistics). Paired-end reads assembly and quality control. Paired-end reads were assigned to individual samples based on unique sample-specific barcodes and truncated by cutting off the barcode and primer sequences. Assembly of paired-end reads was done using FLASH (V1.2.7, https ://ccb.jhu.edu/softw are/FLASH /) 91 . Quality filtering on the assembled sequences (raw tags) was performed to obtain high-quality clean tags setting preset parameters 92 in QIIME (V1.7.0, https ://qiime .org /index.html) 93 . Chimaeras were removed after comparing them with the reference database (i.e., Gold database) using the UCHIME algorithm 94 to detect chimaera sequences, and subsequently, the chimaera sequences were removed, and the useful Tags were finally obtained.
OTU cluster and species annotation. Sequence analysis was achieved using UPARSE software (UPARSE v7.0.1001, https ://drive 5.com/ uparse/) 95 . Sequences with ≥ 97% similarity were allotted to the same OTUs, and the representative sequence for each OTU was screened for species annotation using SSUrRNA database of SILVA Database (https ://www.arb-silva .de/) 96 for species annotation at each taxonomic rank (Threshold: 0.8-1) 97 . To evaluate the phylogenetic relationship of different OTUs, multiple sequence alignment was performed using the MUSCLE software (Version 3.8.31,https ://www.drive 5.com/muscl e/) 98 . Finally, OTU abundance was normalized using the sequence number corresponding to the sample with the least sequences. Downstream analysis of α and β diversity was performed using normalized abundance data.
α-Diversity. The complexity of species diversity (α-diversity) was estimated for each sample using standard indices such as observed-species, sequence depth (Good-coverage) 99 , community richness (Chao1, ACE) 100 , diversity (Shannon, Simpson) 100 . All these indices in beetle samples were calculated with QIIME (Version 1.7.0) 93  β-Diversity. Beta (β) diversity analysis was used to evaluate variances in species complexity in different beetle samples. Beta (β) diversity was calculated in QIIME software (Version 1.7.0) 93 using weighted and unweighted unifrac distances to measure the dissimilarity coefficient between pairwise samples. Non-metric multidimensional scaling analysis (NMDS) 102 , a non-linear model designed for a better representation of non-linear biological data structure, is performed to get principal coordinates and visualize the complex, multidimensional data. Unweighted Pair-group Method with Arithmetic Means (UPGMA) 103 Clustering was performed as a type of hierarchical clustering method to interpret the distance matrix using average linkage and was conducted by QIIME software (Version 1.7.0) 93 . Variation analysis of bacterial community structure between different beetle sample groups was evaluated by standard statistical methodologies such as Analysis of Similarity 104 (ANOSIM, a nonparametric test to evaluate whether variation among groups is significantly larger than the variation within groups), Multi-response permutation procedure (MRPP) 105 analysis (similar with ANOSIM, which aims at determining whether the difference of microbial community structure among groups is significant), ADONIS 106 ( also called permutational MANOVA or nonparametric MANOVA, which is a method of nonparametric multivariate variance test according to the distance matrix, e.g., Bray-Curtis, Euclidean, etc.), AMOVA 107 (similar to ADONIS, which is a kind of nonparametric method aiming at determining whether the difference of microbial community structure among groups is significant). Furthermore, variation analysis of bacterial species between the different beetle sample groups was estimated by T-test 108 and Metastats 109 . Precisely, species with significant intra-group variation were detected via Metastats 109 , rigorous statistical methods based on bacterial abundance measurement. The significance of observed abundances differences among groups is further evaluated via multiple hypothesis-test for sparsely-sampled features and false discovery rate (FDR). Lastly, LEfSe [linear discriminant analysis (LDA) Effect Size] analysis was performed to detect key bacterial species with a significant intra-group variation 110 among tested beetle sample groups. Alternatively, LEfSe is a software aiming at discovering high-dimensional biomarkers and revealing metagenomic features, including genes, metabolic, or taxa, thus can be used to distinguish two or more biological classes. It emphasizes statistical significance, biological consistency, and effect relevance and allows us to identify features of abundance and related classes.
Functional prediction using PICRUSt. Functional prediction of the metagenome was made using PICRUSt 111 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States, version 1.0.0) using OTU table generated in QIIME software (Version 1.7.0) 93 . Kyoto Encyclopaedia of Genes and Genomes (KEGG) 112 were used to calculate the predicted abundance of different gene families. PICRUSt output table was used to build a heat map.