Gut bacterial diversity on the basis of feeding behaviour in different species of thrips (Thysanoptera)

Thrips, tiny to medium sized insects with fringed wings are economically important due to their pestiferous and predatory behaviour. They have a wide range of feeding habitats i.e. leaf feeders, to flower, spore, and fungus feeders. Here, we comprehensively characterized the gut bacteria of nine thrips species collected from two geographical areas of India. Moreover, an attempt to elucidate the impact of different feeding habits (leaf, flower, and spore) on the gut microbiota was also carried out. Comparative analyses revealed that the gut bacterial structure at the phylum level was almost similar with dominant phyla such as Proteobacteria, Firmicutes, Actinobacteria, and Bacteroidetes with few exceptions. Further, at genera level, the gut bacteria were significantly different and revealed the presence of the highest abundance of Rosenbergiella, Wolbachia, Curtobacterium Acinetobacter, and Pseudomonas. Moreover, the gut microbiota was involved in pathways such as Biosynthesis of amino acids, vitamin metabolism etc.


Introduction
Insect gut bacteria play an important role in their development, growth, nutrition, immune system, reproduction, and longevity [1,2]. These microorganisms either protect their host against pathogens, and parasitoids by manufacturing specific toxins or promote the host fitness through several metabolic pathways related to essential amino acids, vitamins, and sterols [3]. Factors such as diet, habitats, and environment play a significant role in influencing the gut bacterial diversity in insects [4]. Earlier studies confirmed this fact that these factors alter the bacterial diversity in Lepidopterans [4][5][6][7][8], Termites [9][10][11], Orthopterans [12], etc. Moreover, researchers also compared the gut bacterial diversity of reared and wild populations and observed significant differences in their gut bacterial diversity [13,14].
Thrips are small, economically important insects under the order Thysanoptera and classified into two suborders Terebrantia and Tubulifera. Presently, 6114 species of thrips are listed across the globe [15], of these 739 are known from India [16]. Thrips are widely distributed in diverse habitats such as leaf, flower, leaf litter, dead branches, decaying bark, fungus, etc. In addition to this, some of the thrips species are pestiferous, gall inducing, predators and Kleptoparasites. The pestiferous species damaged the various agricultural and horticultural crops either by direct feeding or by transmitting plant disease (Topsoviruses) [17]. Keeping in view, the economic importance of this insect group, it is the need of the hour to decipher the gut bacterial diversity along with its predicted functional pathways. Till date, limited studies on gut bacterial diversity of thrips are available in the literature. The gut bacteria of thrips studies started in 2001 with the association of Frankliniella occidentalis and Erwinia sp. [18]. Another study elucidated the mode of gut microbiota transmission in Frankliniella occidentalis and revealed the horizontal transmission of microbiota i.e. egg to adult [19], diet dependent effect of gut bacteria on the host [20], and identification of symbiotic bacteria [21]. Moreover, Dickey et al. [22], assessed the gut bacterial diversity of the most serious pest Scirtothrips dorsalis, and revealed the presence of Actinobacteria and Proteobacteria as the most abundant phyla and Propionibacterium, Stenotrophomonas, and Pseudomonas as the most abundant genera [22].
Recently, Kaczmarczyk et al. [23] reported the gut bacterial diversity in all the developmental stages of the fungivorous thrips (Hoplothrips carpathicus), and revealed the presence of Wolbachia as the most abundant genus in adult and pupa [23]. Subsequently, Gawande et al. [24], reported the microbiome profiling of the onion pest Thrips tabaci from 10 agroclimatic locations in India and observed significant variation in the gut bacterial diversity which may be due to their genetic structure, environment, and habitat [24].
Most of the studies available till date belong either from the lab-reared gut microbiome profiling or are limited to the single species from different environmental habitats as in the case of Thrips tabaci. In our view, the major limitations with the reared specimens were the artificial diet and similar environmental habitat. Thus, the baseline gut microbiota structure of thrips cannot be obtained from the reared specimens as most of the thrips are polyphagous in nature and depends on different vegetation for their feeding. Based on this, in the present work, we made the first largest effort to unveil the gut microbiota structure of nine species of thrips from two geographical locations in India which were further grouped on the basis of their feeding habits such as leaf-feeders, spore-feeders, and flowerfeeders. More specifically the present work will provide the actual baseline data of gut microbiota in thrips for the first time on the basis of environmental and feeding habitats.

Sampling
The thrips specimens were collected from two states of India, West Bengal and Gujarat by beating method onto the white tray and picked with a wetted camel hair brush. During the sample collection, we focused on adults from different feeding habitats from leaves to dead and dry twigs (Table 1). There is no permission required to collect the samples as these are neither endangered nor protected species. Absolute ethanol was used for the preservation of thrips specimens and later stored at 4°C.

DNA isolation, amplification and sequencing
Firstly, thrips specimens were washed individually with PBS solution to reduce the chance of environmental contents. We have given an intersegmental abdominal cut in each specimen to get the gut content in the lysis buffer. After the lysis step, specimens were retrieved and mounted onto the glass slide for morphological identification. Thus, DNA is isolated from 30 specimens individually of each species and pooled for sequencing. Specimens of each species were identified morphologically by Kaomud Tyagi using the available keys [25][26][27][28][29][30][31][32][33]. The lysis solution was further processed for DNA isolation using DNeasy Blood and Tissue Kit (Qiagen) with a manufacturer's protocol. The total DNA of each species was measured using the Qubit 2.0 Fluorometer (Q32866, Thermofisher).

Bioinformatics and statistical analyses
The generated paired end reads of each sample were subjected to bioinformatics analysis. Ninety percent (90%) of the generated paired end reads reaches Q30 and were trimmed 30 bp length for forward reads and 30 bp length for reverse reads and merged based on a minimum overlap of 10 bp to get the 430 bp length. The QIIME2 (ver. 2020.2) was used for merging the forward and reverse reads into a single demultiplexed file [34]. DADA2 pipeline implemented in QIIME2 was used for trimming and de-noising and chimera removal [35]. The non-chimeric reads were assigned into Amplicon Sequence Variants (ASVs). The representative sequence file generated by DADA2 was used for the taxonomic classification using the SILVA database (version 132) with 99% similarity. The online web server Microbiome-Analyst was used for downstream analysis [36,37]. Alpha and Beta diversity were calculated for the normalized data. To measure the diversity richness within the samples, we used two indices (Chao1, Observed) with the statistical method T-test/ANNOVA. Furthermore, to measure both richness and evenness, Shannon and Simpson with statistical method T-test/ANOVA were used. Beta diversity was calculated using Principal Component Analysis (PCoA) and Bray-Curtis index based distance method and PERMANOVA based statistical method. Dendrogram analysis was done for the dissimilarity between samples using the Unweighted Unifrac distance based method and Ward clustering algorithm.
Venn diagram was drawn using a web based tool jvenn to find the unique and shared ASVs between the different feeding habits and environmental habitats [38]. Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) was done for the prediction of the functional metabolic pathways using the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database [39][40][41][42]. The active metabolic pathways between the different feeding habits were plotted using the Stamp software [43]. Nesothrips brevicollis West Bengal Spore-feeding 3.

Bacterial diversity
A total of 378899 raw reads were generated. The sequences obtained from 9 thrips species after the quality trimming, filtering, and de-noising ranged from 10564 (Podothrips sasacola) to 67171 (Hapothrips gowedyi). Further, the rarefaction curve indicated that proper sequencing depth has been achieved in every species for determining the gut microbiota composition. A total of 2157 amplicon sequence variables (ASVs) were obtained, out of which 1286 unique ASVs belong to the species of suborder Tubulifera and 617 unique ASVs belong to the suborder Terebrantia. Both the suborders shared around 254 ASVs towards core gut microbiota. Upon further classification, on the basis of factors i.e. environmental habitat and feeding behaviour, it was observed that thrips species collected from the environmental habitat of West Bengal possess 769 unique ASVs while the species collected from Gujarat possess 1161 unique ASVs (Figure 1a). In addition to this, 227 ASVs resembling core gut microbiota were shared between both environmental habitats. In terms of feeding behaviour, thrips species were classified into three feeding groups i.e. spore-feeders, leaf-feeders, and flower-feeders. Results obtained revealed that they all together shared only 94 ASVs to decipher the core gut microbiota while they possess 505, 732, and 457 unique ASVs (Figure 1b). Further, the highest ASVs were observed in the leaf feeding group, indicating that they have more gut bacterial diversity in comparison to others. Similarly, the gut bacterial diversity of thrips species from Gujarat was more than in West Bengal.

Gut bacterial composition
A total of 11 bacterial phyla were detected with 17 classes, 45 orders, 81 families, and 139 genera in 9 thrips species. The top five abundant phyla in all were Proteobacteria, Actinobacteria, Firmicutes, Bacteriodetes and Cyanobacteria with a mean abundance of 69% (59768 ASVs), 16% (14120 ASVs), 8% (7107 ASVs), 5% (4413 ASVs), 2% (1248 ASVs) respectively. The phylum Proteobacteria was most abundant in the Gujarat environment habitat group due to the flower-feeding species N. gracilipes (95%) and spore-feeding species E. curvipes (76%). Whereas Actinobacteria was highly abundant in West Bengal environment group habitat due to flower-feeding species D. innoxius (54%). Firmicutes were highly abundant in Gujarat environment habitat group and leaf-feeding group due to high abundance in L. furvus (18%). Bacteriodetes were highly abundant in West Bengal environment habitat group and leaf-feeding group due to high abundance in P. sasacola (23%). Cyanobacteria were observed only in West Bengal environment habitat group and sporefeeding group due to high abundance in N. brevicollis (8%) (Figure 1c,d, Table S1).
Total 17 classes were detected which include Gammaproteobacteria (8-93%), Alphaproteobacteria (6-74%), Actinobacteria (1-54%), Bacteroidia (0-21%), and Bacilli (1-19%). The most dominant class is Gammaproteobacteria which was highly abundant in West Bengal environment habitat group and flower-feeding group due to 93% abundance in N. gracilipes, whereas Alphaproteobacteria in Gujarat environment habitat group due to H. gowdeyi with 74% and spore-feeding group in N. brevicollis with 54% abundance. Actinobacteria is dominant in West Bengal environment habitat group due to D. innoxius with 54% and leaf-feeding group due to A. chaetophora with 25%. Bacilli was highly abundant in Gujarat and leaf feeding group due to L. furvus with 19% abundance, Bacteroidia in West Bengal and leaf-feeding group due to P. sasacola with 21% (Table S1).

Alpha and Beta diversity
Alpha diversity based on diversity estimators such as ChaoI, Observed, Shannon and Simpson were carried out. Results obtained based on environmental habitat as a factor, revealed that the thrips species collected from Gujarat region were more diverse than that of the species collected from West Bengal (Figure 2a-d).
Moreover, results obtained based on feeding habits as a factor, revealed that the values of ChaoI and Observed diversity measures for spore-feeders, leaf-feeders, and flower-feeders lie in the range of 120-147, 53-163, 34-145 respectively (Table S2). Based on this, it can be concluded that leaf-feeders were more diverse in comparison to the other feeding groups. On the other hand, for the Shannon and Simpson diversity measures, it was observed that the value for leaf-feeders lies in the range of 3.51-3.  Table S2). Thus, it can be concluded that flower-feeding group has high evenness in its bacterial diversity in comparison to that of the other feeding groups.
The PCoA plot for the Beta diversity analysis based on the Bray-Curtis distance method and PERMANOVA statistical approach using environmental habitat as metadata was presented in the figure (Figure 2i). It was observed that species collected from Gujarat and West Bengal lies in different region with significant variation (p value < 0.05) in the gut microbiota composition. Further, the Beta diversity analysis using feeding behaviour (leaf-feeders, spore-feeders, and flower-feeders) indicated that all thrips species lie in the same region with insignificant (p value < 0.36) variation in gut microbiota composition except N. gracilipes (Figure 2j).
The dendrogram analysis was carried out using the Unweighted Unifrac Distance measure and Ward as a clustering algorithm. The constructed dendrograms based on different feeding habits and environmental habitats indicated similar topologies (Figure 3). The Dendrogram is observed with two clades: clade I with 7 species and clade II with two species (N. gracilipes and P. sasacola). These two species of clade II belonged to West Bengal environment habitat but had different feeding habitats, N. gracilipes from the flower and P. sasacola from the leaf. Clade I is further subdivided into two clades: subclade Ia with 4 species and Ib with 3 species. The subclade Ia species (S. bambusicola, H. gowdeyi, and L. furvus) belonged to Gujarat environment habitat and leaf feeding group. Whereas, clade Ib with 3 species belonged to West Bengal environment habitat and one species A. chaetophora from Gujarat. The clade Ib is further subdivided into Ib 1 with E. curvipes (spore) and D. innoxius (flower); clade Ib 2 into N. brevicollis (spore) and A. chaetophora (flower).

Identification of biomarker using linear discriminant effect size analysis
Linear discriminant effect size (LEfSe) analysis was carried out for the detection of biomarkers associated with different feeding habits, bamboo, leaf, mixed, and spore (Log LDA score 3.0). Ninety-five significant features were detected based on the environmental habitats, wherein, 26 significant features were detected based on the different feeding habits (Figure 4). The phylum Proteobacteria, Actinobacteria, Bacteriodetes, Firmicutes, and Cyanobacteria were identified as potential biomarkers for spore feeding group, Proteobacteria, Choloroflexi, Deinococcus_Thermus for the leaffeeding group. At the genus level, Sphingomonas, Rhizobium, unidentified Cyanobacteria, Microbacterium, Stenotrophomonas, Siphonobacter, Brevundimonas, Streptomonas, Pseudomonas were identified as potential biomarkers for spore feeding group, Acinetobacter, Enhydrobacter, Paracoccus, unidentified Choloroflexi, Deinococcus for leaf feeding group. These taxa identified as biomarkers in the specific feeding group showed high abundance and few are the core gut microbiome.

Prediction of functional metabolic pathways
We investigated the functional metabolic pathways of different feeding habitat groups using PICRUSt2. This analysis inferred the relative abundance of functional genes in 415 metabolic pathways. The carbon metabolism, Biosynthesis of amino acids, Purine metabolism, degradation of aromatic compounds, amino acids metabolism carbohydrate metabolism, cofactors, and vitamin metabolism, were observed in all feeding habitat groups. Out of 415 metabolic pathways 24 pathways were active and shared between flower and spore-feeding group; 3 between leaf and flower; 21 between leaf and spore-feeding group (Figures S1-S3).

Discussion
The gut microbial studies of the order Thysanoptera are very limited and known by only four species (Frankliniella occidentalis, Scirtothrips dorsalis, Thrips tabaci, Hoplothrips carpathicus). Till now, this is the largest effort to explore the gut bacterial diversity and their predicted metabolic pathways from 9 thrips species collected from two geographical locations in India (West Bengal and Gujarat) with three feeding habitats (spore, leaf, and flower). The environmental and feeding habitats of the host species are the most important factors for shaping their gut microbiome, whereas few studies in different groups like Orthopteran revealed that the diet supersede or greatly influenced the gut bacterial diversity [12]. The gut bacterial composition of 9 thrips species revealed 10 bacterial phyla which were represented by 17 classes, 45 orders, 81 families, and 139 genera. Phylum Proteobacteria is most abundant in Gujarat environment habitat species and spore-feeding group species. This result is similar to previous studies in T. tabaci and H. carpathicus, wherein, Proteobacteria was the most abundant phylum, followed by Actinobacteria, Firmicutes, Bacteroidetes, and Cyanobacteria [23]. Proterobacteria is the most dominant in all feeding habitat groups, whereas Actinobacteria is the second most abundant phylum with few exceptions. Cyanobacteria is abundant in one spore-feeding species and two flower-feeding species. This result is similar to an earlier study in T. tabaci with low abundance [24]. The presence of Cyanobacteria also detected in other insects like herbivorous Ensifera and herbivores fishes, indicated that feeding habitats played a greater role in shaping the gut microbiome structure. Cyanobacteria had the capabilities to metabolize carbon and nitrogen [44] and these bacteria were identified in the thrips species of the leaf-feeding group.
At the class level, the Gammaproteobacteria is the most abundant class in all feeding groups followed by Alpaproteobacteria and Bactriodia. Earlier studies revealed that the insect gut is generally governed by the class Gammaproteobacteria which has a high impact on its growth and development [45,46]. Whereas, Alphaproteobacteria was responsible for reproductive regulation, and life history fitness tolerance to the external environment [47].
The gut bacterial structures of thrips species with three feeding habitats and two environmental habitats were significantly different at the genus level. The highly abundant genus Rosenbergiella was detected only in N. gracilipes, but not observed in others feeding groups. This genus was reported in an earlier study of gut bacteria of T. tabaci with 4.5% abundance [24]. Rosenbergiella is commonly known as nectar bacteria and their abundance along with Lactococcus, and Asaia had been studied to evaluate the effect on the longevity of aphid parasitoids [48]. It revealed that these bacteria do not affect longevity but have a strong impact on the foraging behaviour of animals along with plant fitness. Pseudomonas was the most abundant genus in spore feeding groups and was also reported in earlier thrips [22,24] and other insect gut microbiome studies like Lepidoptera [49].
Stenotrophomonas (cellulolytic bacteria) was the most abundant genus in P. sasacola and S. bambusicola of the leaf-feeding group which usually inhabits on bamboo. Earlier, this genus was reported in the wood-boring beetles, Saperda vestita and Agrilus planipennis, and bark beetle Dendroctonus rhizophagus, wherein these were involved in the degradation of cellulose [50][51][52].
Curtobacterium known as latex degrading bacteria was the most abundant in flower-feeding habitats, this result is similar to an earlier report on lepidopteran species, Hyles euphorbiae feed on Euphorbia sp. and Brithys crini on Pancratium maritimum [53]. The genus Acinetobacter is majorly reported in all feeding habitat groups and also reported in Lepidoptera, Lymantria xylina is known as a major defoliator of hardwood and fruit trees [54]. Stenotrophomonas, Pseudomonas, and Acinetobacter were detected in all feeding habits, involved in the fermentation, degradation, and oxidation of aromatic compounds [52]. Pseudomonas and Acinetobacter also play a significant role in the nutrient supplements in sawflies [55]. The genus Methylobacterium was detected in all feeding habitats except fungus as reported in other lepidopteran insects and potentially contributed towards nitrogen fixation [56].
Wolbachia is known as the most prevalent endosymbiont, which alters the sex ratio and induces cytoplasmic incompatibility and parthenogenesis in 40% of arthropod species [57][58][59]. Remarkably, Wolbachia was the dominant genus in a leaf-feeding group with high abundance in H. gowdeyi (62%), P. sasacola (3%), and low in flower-feeding species N. gracilipes (1%). This genus was not observed in other feeding groups in contrast to earlier leaf-feeding species like T. palmi [60] and T. tabaci [24]. During the collection of these species which were detected with Wolbachia in the gut, we were not able to collect a single male specimen. This result directly indicated that Wolbachia induces cytoplasmic compatibility or male killing or parthenogenesis. Therefore, it can be concluded that Wolbachia might be used as a controlling agent for thrips pest species. Alpha diversity indicated that the diversity and richness of the gut bacteria appeared to be higher in Gujarat than in West Bengal. In respect to the feeding habits, the leaf-feeding group has the highest gut bacterial diversity followed by spore, and flower indicating that the diversity and richness of gut bacteria in insects was higher which usually feed on leaves of various plant families. Dendrogram analysis revealed that the separate clade of P. sasacola + N. gracilipes from the rest of the species, in contrast to their feeding habitat and environmental habitat was due to the presence of Rosenbergiella in N. gracilipes, Candidatus, and Alkanindiges in P. sasacola. PCoA analysis also indicated the same results for N. gracilipes and P.sasacola, as these two species were clustered separately from the rest. The cladding of S. bambusicola, H. gowdeyi, and L. furvus occurred due to the presence of the Wolbachia genus. Thrips species with different environmental habitats showed differentiation in the gut bacterial diversity. In contrast, the gut with their feeding habitats did not show a clear differentiation except for N. gracilipes. This result indicated that not only feeding behaviour but the environment was also playing a significant role in shaping the gut microbiome structure.

Conclusion
The present work is the first comprehensive study unveiling gut bacterial diversity based on two parameters (environmental and feeding habitats) of nine thrips species. The result on both parameters at the phylum level was almost similar in all but varied at the genus level. Moreover, the richness of gut bacteria in the leaf feeding groups was much higher in others, whereas evenness was higher in flowers than in others. Further, unique bacteria for the flower-feeding, and the leaffeeding groups can be used as a potential biomarker. There is no unique genus for the spore feeding group. The gut bacteria are involved in the Biosynthesis of amino acids, degradation of aromatic compounds, carbon fixation, Methane metabolism, and other cofactors and vitamin metabolism. This study indicated that the different feeding habits and environmental factors may have a great effect on shaping the gut bacterial structure. Furthermore, in order to understand the pattern of gut bacterial diversity in different feeding habitats and the relationship with the environment in this group, more sampling and study of the gut bacteria is essential.

Authors contribution
KT, IT, and VK were involved in Conceptualization, Data Curation, and Methodology; KT, DS, and AP were involved in specimen collection, identification, and DNA isolation; KT, and IT were involved in Software and Bioinformatics analysis; IT, VK, DS and AP involved in functional analysis; All the authors were involved in Manuscript writing and Editing. VK supervised the project and Funding acquisition.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by ZSI Core funding, Ministry of Environment Forest and Climate Change, New Delhi, Republic of India.

Availability of data and material
The generated raw reads were submitted under the Bio Project ID RJNA782448 to the National Centre for Biotechnology Information (NCBI) GenBank Portal.