Does diet breadth affect the complexity of the phytophagous insect microbiota? The case study of Chrysomelidae.

Chrysomelidae is a family of phytophagous insects with a highly variable degree of trophic specialization. The aim of this study is to test whether species feeding on different plants (generalists) harbour more complex microbiotas than those feeding on a few or a single plant species (specialists). The microbiota of representative leaf beetle species was characterized with a metabarcoding approach targeting V1-V2 and V4 regions of the bacterial 16S rRNA. Almost all the analysed species harbour at least one reproductive manipulator bacteria (e.g., Wolbachia, Rickettsia). Two putative primary symbionts, previously isolated only from a single species (Bromius obscurus), have been detected in two species of the same subfamily, suggesting a widespread symbiosis in Eumolpinae. Surprisingly, the well-known aphid symbiont Buchnera is well represented in the microbiota of Orsodacne humeralis. Moreover, in this study, using Hill numbers to dissect the components of the microbiota diversity (abundant and rare bacteria), it has been demonstrated that generalist insect species harbour a more diversified microbiota than specialists. The higher microbiota diversity associated with a wider host-plant spectrum could be seen as an adaptive trait, conferring new metabolic potential useful to expand the diet breath, or as a result of environmental stochastic acquisition conveyed by diet.

Introduction

Insects are colonized by a variety of microorganisms, prevalently living as commensals, but which in many cases can confer either beneficial or detrimental effects to their host (e.g., Douglas, 2009; Kikuchi et al., 2012; Engel and Moran, 2013; Clay, 2014; Douglas, 2015; Hurst and Frost, 2015; Wang et al., 2020). Since symbiont‐mediated traits highly influence the host nutrition, in herbivorous insects this influence is often crucial in the interaction with the host plant (Hansen and Moran, 2014; Giron et al., 2017; Mason et al., 2019; Frago et al., 2020; Mason, 2020). Most of the microorganisms that can be found within the insect body colonize the gut lumen, but they can be also hosted in specialized organs often connected to female genitalia, especially when the vertical transmission of the symbiont is required (Stammer, 1935, 1936; Mann and Crowson, 1983; Becker, 1994). The most specialized bacterial symbionts live inside the insect cells are vertically transmitted, and show drastic genome reduction (e.g., Blattabacterium, Buchnera), usually maintaining only the metabolic pathways involved in providing functional traits to the host (Boscaro et al., 2017; Latorre and Manzano‐Marín, 2017; Ankrah et al., 2018). Other bacteria are able to colonize insect cells, such as the reproductive manipulators belonging to the so‐called male‐killing group (e.g., bacteria of the genera Wolbachia and Rickettsia) that can manipulate the host reproduction to maintain their infection across generations and spread within the population (Harris et al., 2010; Correa and Ballard, 2016; Larracuente and Meller, 2016).

Most studies investigating the relationship between bacterial symbionts and the insect host have been conducted on model species, mainly focusing on single interactions. More recently the advent of next‐generation sequencing techniques coupled with the 16S rRNA‐based approach for the bacterial taxonomy has greatly facilitated the characterization of the full microbiota associated with non‐model organisms, allowing to expand the experimental scale (e.g., Montagna et al., 2015a; Mohammed et al., 2018; Ziganshina et al., 2018; Kolasa et al., 2019). This innovation opened the possibility to characterize the microbiota associated with several wild species and so to investigate the correlations between the composition of microbial communities and several ecological or physiological traits of the insect host, such as the breadth of the insect diet (Colman et al., 2012; Yun et al., 2014). Indeed, insects feeding on several plant species may be expected to harbour more complex microbial communities. The gut microbiota composition can be influenced by the diet, directly since food may inoculate bacteria able to colonize the insect gut or indirectly by promoting the growth of specific bacteria (Pérez‐Cobas et al., 2015; Montagna et al., 2015b; Chouaia et al., 2019; Muturi et al., 2019). So, insects with a wider food source spectrum are expected to be colonized by a higher diversity of microbial taxa. Anyway, the higher diversity in the microbiota of generalist species could also be due to the wider metabolic potential, conferred by the presence of a more variegated microbial community, which makes those insects able to exploit several different food sources. The covariation of microbiota diversity and breadth of the animal diet has been investigated also in non‐insect taxa, usually achieving inconclusive results that do not support the hypothesis of a higher diversity in the microbiota of generalist species (e.g., Kartzinel et al., 2019; Chen et al., 2021).

Leaf beetles (Coleoptera: Chrysomelidae), including ~40 000 species worldwide, constitute one of the most diverse insect groups in the world (Leschen and Beutel, 2014). This Coleoptera family includes almost only phytophagous species, feeding on leaves or other plant organs at least at the adult stage. The degree of trophic specialization is highly variable, since some leaf beetle species can exploit only one or few specific plant species as a food source, while others can feed on hundreds of plant species belonging to several different families. This makes Chrysomelidae a perfect model to investigate the relationship between the level of microbiota complexity and the breadth of the host plant spectrum. Furthermore, leaf beetles are of great interest for the presence of vertically transmitted symbionts (i.e., in Donacinae, in Cassidinae and in Eumolpinae), which are harboured in specialized host organs associated with gut and genitalia (Stammer, 1935, 1936; Tayade et al., 1975; Mann and Crowson, 1983; Becker, 1994). Bacteria of the genus ‘Candidatus Macropleicola’ (Enterobacteriaceae) are hosted in specialized organs at the midgut–hindgut junction of Donacinae. These bacteria show a tight co‐speciation with the insect host and are involved in supporting its nutrition providing essential nutrients during the larval stage (essential amino acids, riboflavin) and digestive enzymes (pectinases) to the adult insects (Kölsch et al., 2009; Kölsch and Pedersen, 2010; Kleinschmidt and Kölsch, 2011; Reis et al., 2020). Similarly, ‘Candidatus Stammera capleta’ (Enterobacteriaceae), hosted in specialized organs associated with the foregut of several Cassidinae species, is involved in pectinase production (Salem et al., 2017, 2020). Within Eumolpinae only a single species is known to host symbionts in specialized organs, Bromius obscurus (Stammer, 1936). This symbiosis has been less studied, but two different symbionts have been isolated from it (Kölsch and Synefiaridou, 2012). The first one (henceforth B. obscurus symbiont A) is hosted intracellularly in blind sacs at the foregut–midgut junction and extracellularly in female‐specific genital accessory organs, suggesting the presence of vertical transmission. The second one (henceforth B. obscurus symbiont B) is hosted in small crypts at the end of the midgut and is phylogenetically related to bacteria species living in the gut lumen, not tightly associated with the host (Kölsch and Synefiaridou, 2012; Fukumori et al., 2017).

Previous studies on the bacteria associated to leaf beetles were mainly focused on reproductive manipulators (Clark et al., 2001; Keller et al., 2004; Kondo et al., 2011; Roehrdanz and Wichmann, 2013; Montagna et al., 2014; Krawczyk et al., 2015; Kolasa et al., 2017; Takano et al., 2017; Gómez‐Zurita, 2019), single species of economic importance (Muratoglu et al., 2011; Chung et al., 2013; Ali et al., 2019; Ludwick et al., 2019; Wang et al., 2019; Shukla and Beran, 2020) or few strictly related species (Kelley and Dobler, 2011; Montagna et al., 2015a; Blankenchip et al., 2018; Wei et al., 2020). The present study aims to characterize the microbiota associated with a selection of leaf beetle species, representative of the taxonomic diversity and of the various degrees of trophic specialization. In detail, it aims: (i) to determine the principal bacterial taxa that characterize the microbiota of the selected leaf beetle species, also detecting the presence of important insect symbionts (e.g., Wolbachia) and symbionts typically present in specific Chrysomelidae subfamilies (e.g., Donacinae, Cassidinae); (ii) to test the hypothesis that the microbiota of generalist phytophagous species is more complex than the microbiota of more specialist species.

Results

Efficiency and taxonomic resolution of the 16S rRNA gene V1–V2 and V4 regions

From the 30 Chrysomelidae species analysed a total of 841 822 (mean per sample = 28 060.7) and 1 711 075 (mean per sample = 57 035.8) raw reads have been obtained from the sequencing of the V1–V2 and V4 regions of the bacterial 16S rRNA respectively. Raw sequences have been deposited on the NCBI SRA database under the project accession number PRJNA729224. After the denoising and filtering steps, the V1–V2 dataset consisted of 1080 amplicon sequence variants (ASVs) (total reads = 449 368; mean per sample = 14 978.9) and the V4 dataset consisted of 1572 ASVs (total reads = 1 047 226; mean per sample = 34 907.5). All the ASVs assigned to mitochondria (<0.01% of the V1–V2 reads, 0.74% of the V4 reads) or chloroplast (45.3% of the V1–V2 reads, 17.8% of the V4 reads) have been excluded from further analyses. Regarding the taxonomic identification of the ASVs, 119 bacterial genera have been identified by both regions while 162 genera are present only in the V4 dataset and 35 genera only in the V1–V2 dataset. Comparing the results of the taxonomic assignment of the two regions, the V4 region results the marker almost always more efficient in detecting bacterial taxa (Supplementary Fig. 1), with few exceptions in which those for the V1–V2 region slightly outperform the others (e.g., the genus Brevundimonas and Aeromonas). The estimated diversity using the two regions separately (Supplementary Fig. 2) is identical when putting much weight on the most abundant species (q = 2), while the V4 region provides slightly higher estimates when increasing the weight of rare species (q = 1, q = 0).

Microbiota composition

The most represented bacterial classes associated to the analysed Chrysomelidae species (Figs 1A and 2) are Alphaproteobacteria (~39%), Gammaproteobacteria (~45%) and Bacilli (~14%). Several genera belonging to Bacteroidia are also present in the microbiota of the selected species (Fig. 2) but this class constitutes only ~1% of the total dataset. Within Alphaproteobacteria the most abundant genera recorded are Wolbachia, Rickettsia and Sphingomonas. Wolbachia is the most represented genus in the dataset (~30% of the total reads) and sequences belonging to this genus have been found in all the species except Chrysomela saliceti and Timarcha tenebricosa (Fig. 1B, Supplementary Table 1). In most cases Wolbachia sequences represent a low percentage of the sample reads (<1%), while in eight species it comprises the most abundant ASVs. Rickettsia is another well‐represented genus (Fig. 1B, Supplementary Table 1). Reads assigned to Rickettsia have been found in nine species and represents a high percentage of the reads from Hispa atra (~31%) and from all the sampled species belonging to Clytrini tribe (Labidostomis longimana ~18%, Clytra quadripunctata ~44%, Smaragdina affinis ~92%). Sphingomonas is also quite common in the dataset (Fig. 1B). It is present in all the sampled species, with the only exception of Plateumaris consimilis, and it reaches the highest densities in Timarcha tenebricosa (~10%) and Cryptocephalus transcaucasicus (~30%). Within Gammaproteobacteria the most abundant genus is Pseudomonas (Fig. 1A). Pseudomonas is the second most abundant genus in the dataset (~12% of the total reads) and it is the only bacterial genus that is present in all the sampled species, with relative abundances that varies from 0.3% to 82%. In total 141 ASVs (corresponding to 15 97%‐similarity OTUs) are assigned to Pseudomonas, 69 of them are present only in one species so that 21 species have at least one unique ASV assigned to this genus. Another surprisingly well‐represented bacterial genus belonging to Gammaproteobacteria is Buchnera (Fig. 1A). It has been found in 21 species, mostly at low abundance but representing a quite high proportion of the reads obtained from Orsodacne humeralis (~22%). The most abundant genus belonging to Bacilli is Spiroplasma, followed by Paenibacillus and Brevibacillus (Fig. 1A). Spiroplasma is the dominant genus in Crioceris paracentesis (~98%) and it has been found also in 18 other species, but with low densities (Fig. 1B). While in the microbiota of Lilioceris merdigera the dominant genus is Paenibacillus (~58%), that is also present in few other species but with low abundances.

Fig. 1

Taxonomic composition of the bacterial microbiota of Chrysomelidae.

A. Bar plot representing the composition of the microbiota of each Chrysomelidae species at the genus level (or any higher taxonomic level when the identification at the genus level was not possible). Colours represent different bacterial ranks, as reported in the legend, and the height of each box corresponds to the average relative abundance of each bacterial rank. Only bacterial ranks representing on average at least 5% of the reads in one species are shown, less abundant bacteria are included in the group ‘Others’.

B. Relative abundance of reproductive manipulators (Wolbachia, Rickettsia, Spiroplasma) and bacterial symbionts present only in Chrysomelidae. Colours representing different bacterial genera are reported in the legend (in bold). Only insect species with a relative abundance of these bacteria of at least 3% are shown.

Fig. 2

Diversity and composition of the Chrysomelidae microbiota. On the left side the ancestral state reconstruction of the microbiota diversity (Hill numbers, q = 1) is plotted as a colour gradient on the ML phylogenetic tree of the selected Chrysomelidae species. The heatmap represents the relative abundance of each genus with colours corresponding to the bacterial class (only classes with at least five different bacterial genera are reported, other classes are included in the category ‘Other bacterial classes’). On the right side, together with the insect species names, the presence of one or more reproductive manipulators or Chrysomelidae specific symbionts (one or two bacterium icons) and the trophic classification of the insect (specialists with one leaf icon; generalists with two leaf icons) are reported.

The results of the NCBI blast search (Supplementary Table 2) and phylogenetic tree inference (Supplementary Fig. 3) shed further light on the bacterial taxa that characterize the microbiota of Chrysomelidae. Two ASVs, which have been found only in Cryptocephalus fulvus (53.6%), have been assigned by the naïve Bayes classifier (confidence >0.95) only to the domain level and the top hits of the blast search (query coverage 100%, identity >80%) correspond to uncultured bacteria isolated from acidic biofilm in caves (DQ499258). A huge number of sequences obtained from Donacinae (Donacia obscura 55.3%, Plateumaris consimilis 93.8%) belong, with high confidence (CP046230; query coverage 100%, identity 99.7%), to the vertically transmitted endosymbiont widespread in this subfamily (Kölsch et al., 2009). In both the maximum likelihood (ML) trees those sequences cluster with sequences obtained from the bacterial symbiont of Donacinae with quite high confidence (bootstrap values >70). Similarly, most of the sequences obtained from Cassidinae (Cassida inopinata 74.7%, Hypocassida subferruginea 93.1%) resulted to belong to ‘Candidatus Stammera capleta’ (CP024013; query coverage 100%, identity 98.8%) and cluster with sequences of this species in the ML trees with high confidence (bootstrap value =100). Three ASVs from the V1–V2 region assigned to Enterobacterales are present only in Macrocoma henoni (23.9% of the reads). Blast search top hits (query coverage 100%, identity >78%) correspond to endosymbionts of weevils (AP018159, KX067892) while in the ML tree they cluster with a sequence from the B. obscurus symbiont A (LC273302) with high confidence (bootstrap value =99). Also, two ASVs from the V4 region from M. henoni (that represent almost all the reads previously assigned to Buchnera in this species) clustered together with sequences of the B. obscurus symbiont A (bootstrap value =88; Supplementary Fig. 3); the blast search confirms this taxonomic annotation (query coverage 100%, identity 93.9% and 93.5% with LC273302 and JQ805030 respectively). Six ASVs are present only in Chrysochus asclepiadeus and represent 93.3% of the reads from this species, one of them can be assigned to the B. obscurus symbiont A (LC273302, JQ805030; query coverage 99%, identity >87%, bootstrap value =97). Among the remaining ASVs blast search top hits assigned two ASVs to Lelliottia amnigena (LR134135; query coverage 100%, identity >98%) and the other three ASVs to Klebsiella sp. (MN860163, LR134475, MT279983, MT255043; query coverage 100%, identity >98%). In the ML trees (Supplementary Fig. 3) those sequences are part of a clade that includes Klebsiella and Lelliottia, but also other bacterial genera including symbionts of Hemiptera (JQ322760, HM156667, AB650515, AY620432) and the B. obscurus symbiont B (JQ805033).

Microbiota diversity

Diversity estimates for the group of generalist species, identified as those feeding on several plant families, are always higher than estimates for specialist species (Fig. 3; Supplementary Fig. 2). Similar results have been obtained also defining the two trophic groups (i.e., generalists, specialists) by working at the level of plant genera (Supplementary Fig. 4). As an example, with q = 2 (i.e., counting mainly the dominant taxa) the diversity estimated for generalist species is almost twice the diversity estimated for specialist species (coverage >0.3). Also, the diversity partitioning analysis in the framework of Hill numbers confirms this pattern (Table 2). The α‐diversity component (average diversity of single species microbiotas) is always higher in generalist species, regardless of the value of the order parameter (q). Also, the γ‐diversity (diversity of the microbiota of all the species together) is higher in generalist species, except in the case of q = 0 (i.e., counting mainly the rare species). When using an intermediate weight (q = 1; counting mainly the common species) the gamma diversity estimated for the entire dataset (5.5) has an intermediate value between specialist species (4.9) and generalist species (6.4), as expected, while using other values for the order parameter there are no clear differences. The β‐diversity estimates (γ‐diversity/α‐diversity; corresponding to differences between samples) are similar in generalists and specialists, except for q = 0.

Fig. 3

Microbiota diversity of specialist and generalist Chrysomelidae defined using the plant taxonomic level of family (specialists feed on plants all belonging to the same family, generalists feed on plants belonging to different families). Coverage based rarefaction/extrapolation curves of the Hill numbers estimated for three values of the order parameter (q = 0, q = 1, q = 2). The x‐axis represents the coverage (that estimates the completeness of the sampling) and the y‐axis represents the Hill number estimates, 95% confidence interval is also reported. As reported in the legend, colours correspond to the trophic category (specialist or generalist) and line type to the methodological approach (interpolation or extrapolation).

Table 2

Diversity partitioning.

Group	q	α‐diversity	γ‐diversity	β‐diversity
Total dataset	0	35.2	318.4	9.0
	1	1.9	5.5	2.9
	2	1.2	3.0	2.5
Specialists	0	26.8	318.4	11.9
	1	2.7	4.9	1.8
	2	2.0	2.9	1.4
Generalists	0	39.1	318.4	8.1
	1	3.5	6.4	1.8
	2	2.3	3.0	1.3

The three diversity components (α‐diversity, within sample diversity; γ‐diversity, whole group diversity; β‐diversity, among samples diversity) are reported for the whole dataset (total dataset) and for each of the two trophic categories considered (specialists and generalists).

The ancestral state reconstructions of the microbiota diversity estimates along the Chrysomelidae phylogenetic tree show no clear pattern (Fig. 2), in fact, no phylogenetic signal has been recorded. Some phylogenetic clades share low levels of microbial diversity (e.g., Donacinae, Cassidinae) but it is probably related to the presence of primary symbionts at high abundances (Fig. 2). In fact, most of the species hosting primary symbionts (e.g., ‘Candidatus Stammera capleta’) or reproductive manipulators (e.g., Wolbachia) have low microbial diversity estimates (especially for q < 2). In any case, there are still several species hosting reproductive manipulators that have highly diverse microbiota (e.g., Labidostomis longimana, Chaetocnema hortensis).

Discussion

The V4 region of the 16S rRNA, one of the two markers used in this study to characterize the microbiota associated to Chrysomelidae, allowed to obtain a higher number of ASVs and seems less prone to chloroplast contamination in respect to the second marker adopted, the V1–V2 region. Moreover, several bacterial genera were identified only through the examination of V4 region reads (Supplementary Fig. 1) and the diversity analyses suggest that this region better recovers rare taxa (Supplementary Fig. 2). These results are in accordance with what was found in previous studies supporting the use of the V4 region of the 16S rRNA gene in metabarcoding studies on bacteria (Zhang et al., 2018; Chen et al., 2019). Nevertheless, the V4 region missed the amplification of some bacterial taxa (e.g., Brevundimonas and Aeromonas) and the values of the diversity indices obtained combining the two regions resulted higher than those obtained from a single region (Supplementary Fig. 2). These results support the use of multiple marker regions to increase the resolution of metabarcoding studies targeting bacterial communities.

Within the microbiota associated with the 30 species of Chrysomelidae analysed in this study, three endosymbiotic bacterial genera belonging to the so‐called male‐killing group (Engelstädter and Hurst, 2009) have been identified: Wolbachia (in 28 species), Spiroplasma (in 17 species) and Rickettsia (in nine species) (Supplementary Table 1). For the majority of the analysed Chrysomelidae, the association with these bacteria is reported in this study for the first time. The presence of reproductive manipulators, such as Wolbachia, in Chrysomelidae is well known (Montagna et al., 2014; Kajtoch and Kotásková, 2018; Gómez‐Zurita, 2019). These endosymbiotic bacteria are usually abundant in infected species, tending to dominate the community, as observed in this study for four Chrysomelidae species, where Wolbachia represent 94%–99% of the reads. Endosymbionts can also represent only a minimum fraction of the bacterial community, e.g., Wolbachia represent less than 0.05% of the reads in nine species analysed in this study. The latter cases could be signs of horizontal acquisition that did not lead to an infection (Rasgon et al., 2006; Pietri et al., 2016; Chrostek et al., 2017; Kolasa et al., 2017; Cardoso and Gómez‐Zurita, 2020) rather than real infections able to produce effects on the host even with a low bacterial titre (Richardson et al., 2019). Surprisingly, Orsodacne humeralis was found to host Buchnera (Gammaproteobacteria: Enterobacteraceae) representing the 22.1% of bacterial reads obtained for this species. This bacterium is a well‐known endosymbiont that is strictly associated with aphids (Buchner, 1965; Shigenobu and Wilson, 2011) and to our knowledge infections caused by it in a non‐aphid host have never been reported. For this reason, it is also hard to determine the relationship between Buchnera and O. humeralis microbiota (e.g., acquisition from the environment, commensality, symbiosis). Pseudomonas is the second most abundant bacterial genus found to be associated with the Chrysomelidae species of this study and it is the only one detected in all analysed species. Species of this genus can live under diverse environmental conditions; they are ubiquitous in soil, water and are important pathogens of plants and animals (Moore et al., 2006). Pseudomonas species are also known to play a functional role in insect symbiosis (e.g., providing digestive enzymes) (Piel et al., 2004; Huang et al., 2012; Ceja‐Navarro et al., 2015; Zhang et al., 2020a,b). The high variety of Pseudomonas species makes it difficult to distinguish among environmental contamination (presumably from the food source), facultative association and functional symbiosis. Nevertheless, the high prevalence and uniqueness of Pseudomonas ASVs in some samples allows to suppose that, at least in these cases, it could represent a symbiont potentially playing a functional role for some Chrysomeldiae species. Three Chrysomelidae subfamilies (Donacinae, Cassidinae, Eumolpinae) are known to host specific bacterial symbionts in specialized organs associated with the gut. Symbionts of Donacinae and Cassidinae have been intensively studied in the last years and are known to support host nutrition supplying digestive enzymes and/or providing essential nutrients lacking in the insect diet (Kleinschmidt and Kölsch, 2011; Salem et al., 2017, 2020; Reis et al., 2020). While for Eumoplinae those kinds of symbiosis have been less studied and are known only in B. obscurus. In both the species of Donacinae included in this study the microbiota is dominated by an endosymbiont already known to be widespread within the species of the subfamily (Fig. 1B, Supplementary Table 1). Specifically, in Plateumaris consimilis ~94% of the sequences are assigned to the symbiont isolated from that same species in Reis et al. (2020), while in Donacia obscura ~55% of the sequences have been assigned to the symbiont isolated in the same study from Donacia cinerea and Donacia marginata, since no reference sequences are available for Donacia obscura symbiont. Most of the reads obtained from both the species of Cassidinae included in this study, Cassida inopinata (74.7%) and Hypocassida subferruginea (93.1%), have been assigned to symbiont of Cassidinae ‘Candidatus Stammera capleta’ (Stammer, 1936; Salem et al., 2017) (Fig. 1B, Supplementary Table 1). This result is the first report of ‘Candidatus Stammera capleta’ in these two Cassidinae species. The two analysed species of Eumolpinae (M. henoni and Chrysochus asclepiadeus) host bacterial taxa previously reported only from B. obscurus (B. obscurus symbiont A and B). The microbiota of M. henoni is dominated by the B. obscurus symbiont A. Since in B. obscurus this intracellular bacterium is present in specialized gut organs and in the female genitalia, it is possible to hypothesize a similar localization and a vertical transmission mechanism also in M. henoni but further studies are needed to investigate this symbiotic relationship and confirm this hypothesis. The symbiont A of B. obscurus is also present in Chrysochus asclepiadeus, but with a low abundance (~2%). The microbiota of Chrysochus asclepiadeus is dominated by a group of closely related bacterial taxa including the B. obscurus symbiont B together with the two bacterial genera Lelliottia and Klebsiella. Lelliottia spp. are usually isolated from plants, water and clinical samples (Brady et al., 2013). This genus has been also found in insect microbiota but usually at low abundances (e.g., Wang et al., 2019; Xu et al., 2019). Klebsiella spp. are often isolated from a variety of environmental sources such as soil, vegetation, water and animals (Brisse et al., 2006), but some species are also known to play functional roles in insect symbiosis (e.g., providing enzymes and antibiotics) (Dillon et al., 2002; Dantur et al., 2015; Miyashita et al., 2015). The high prevalence of ASVs closely related to the B. obscurus symbiont B in Chrysochus asclepiadeus suggests the presence of similar symbioses in closely related Eumolpinae species, which should be further investigated.

Interestingly, in the 16S rRNA phylogenetic analyses (Supplementary Fig. 3) ‘Candidatus Stammera capleta’, the endosymbiont of Donacinae and the B. obscurus symbiont A (together with five M. henoni ASVs and one Chrysochus asclepiadeus ASV) clustered in a clade with the most specialized Enterobacteriaceae symbionts (e.g., Buchnera, Blochmannia, Nasonia, Baumannia). However, the symbiont B of B. obscurus is placed in a separate clade with more generalist bacteria (e.g., Escherichia, Enterobacter, Serratia, Klebsiella, Lelliottia). This supports the hypothesis that the B. obscurus symbiont B, also present in Chrysochus asclepiadeus, is related to quite generalist gut bacteria, suggesting a loose association with the host. While the B. obscurus symbiont A is in the same clade with the other subfamily specific Chrysomelidae symbionts, and since it has been detected also in both the Eumolpinae species included in this study, it is probably widespread in the subfamily.

The comparison of α‐diversity metrics between the specialist and the generalist species included in this study confirms that those usually feeding on several plant families harbour a more diversified microbiota (Fig. 3). The higher microbiota richness of generalist insects has been previously observed comparing different insect orders and broad diet categories (e.g., detritivores vs. herbivores/carnivores) (Colman et al., 2012; Yun et al., 2014). Also, a study performed on Tephritidae supports this hypothesis (Ventura et al., 2018), while other studies on different taxonomic groups do not confirm this pattern (Blankenchip et al., 2018; Rothman et al., 2020). The higher microbiota diversity observed in generalist insects can be the result of bacteria randomly acquired from the environment, without any specific functional role in the host physiology, or due to the establishment of a diversified microbiota that provides adaptive advantages to the host (i.e., a wider metabolic potential that allows the exploitation of diversified food sources). An exemplar case can be identified in detritivorous insects that harbour one of the richest microbiotas among insects (Colman et al., 2012). Detritivorous insects' food source is composed of substrates of different origins that are colonized by a high variety of bacterial taxa, potentially contributing to insect microbiota diversity. On the other hand, detritus includes some of the most difficult molecules to digest (e.g., lignocellulose), thus insects could benefit from the amplified metabolic capabilities supplied by a richer microbiota. Also in the case of phytophagous insects, the higher microbiota richness observed in generalist species can be easily related to the acquisition of different bacteria that are part of the environmental microbial communities (Montagna et al., 2015b; Chouaia et al., 2019; Hannula et al., 2019; Jones et al., 2019). The host plant–soil system is one of the major drivers of the phytophagous insect microbiota (Hannula et al., 2019), thus feeding on more than one plant species can highly influence insect's microbiota diversity and composition (Jones et al., 2019). Secondary metabolites (Zhang et al., 2020a; Zhang et al., 2020b) and plant defences (Chung et al., 2017) often play a fundamental role in these plant–insect–microbiota interactions. Moreover, bacteria colonizing plant surfaces and tissues can be able to degrade the toxic compounds produced by the plant itself (e.g., Shukla and Beran, 2020; Leite‐Mondin et al., 2021). So, insects feeding on plants can acquire bacteria that, if established in their microbiota, can provide adaptative advantages. In fact, a richer microbiota determines a wider range of metabolic capabilities that can help the phytophagous insect to overcome the defences of different host plants (e.g., Martinez et al., 2019; Santos‐Garcia et al., 2020), thus also allowing the expansion of the trophic spectrum. Distinguishing between the processes that shape the microbiota diversity of generalist insects, especially if phytophagous, is quite difficult. Indeed, our results suggest that probably both the random acquisition from the environment and the adaptive advantage of an amplified metabolic potential participate to increase the diversity of the microbiota of generalist species. In this study, the estimated diversity of the microbiota is always higher in generalist species. This is observed both when the diversity value is estimated considering all bacteria (the weight is mostly on rare and low‐abundance species, more likely acquired from the environment; q = 0), as well as when the weight is on more common species (considering mainly medium‐high abundance bacteria having a possible functional role; q = 1, q = 2) (Fig. 3). The importance of the food source in influencing the composition of the microbiota is also highlighted by the higher β‐diversity observed in specialist insects when focusing on rare bacteria (q = 0) (Table 2). In fact, the microbiota of each specialist species results simpler than that of each generalist species (lower α‐diversity) but considering together all the species in each of the two groups the overall diversity reaches similar levels (same γ‐diversity). Based on previous results, the overall microbiota diversity of the specialists is mainly due to the high amount of not‐shared bacterial taxa among species (reflected by a high β‐diversity). This can be explained by the acquisition of phylogenetically distant bacteria from the host plant exploited by each insect species, supporting the importance of the food source in influencing the microbiota of phytophagous insects. The higher α‐diversity of the microbiota harboured by generalist species (approximately twice higher than that of specialists) is confirmed also when focusing on dominant bacteria (q = 2; Fig. 3). In this last case the difference is probably related to non‐transient bacteria that may have a functional role in insect physiology. These results support the hypothesis that the high bacterial diversity hosted by generalist insects can expand the host metabolic potential enabling the exploitation of different food sources. Further studies, with an increased sample size or focusing on other phytophagous insect groups, are needed to confirm the pattern here observed and to better clarify the most influential causes of this phenomenon.

Experimental procedures

Species selection and host plant information

Thirty adult insects, collected from vegetation by sweep net and identified as belonging to 30 different species of Chrysomelidae, have been selected for this study (Table 1). The selection was performed to maximize the taxonomic coverage and the representativeness of the variability in the trophic specialization. The sampling includes representatives of the 10 main subfamilies of Chrysomelidae: Alticinae, Chrysomelinae, Galerucinae, Donacinae, Criocerinae, Cassidinae (including Hispini), Cryptocephalinae (including Clitrini), Eumolpinae, Orsodacninae, Zeugophorinae. The list of plants included in the diet of each Chrysomelidae species was compiled from a database of the host plants of Euro‐Mediterranean Chrysomelidae (Magoga et al. in preparation). The trophic spectrum of the selected species ranges from exclusively monophagous species, restricted to feed on a single plant species (e.g., Chrysocus asclepiadeus feeds only on Vincetoxicum hirundinaria, Apocynaceae), to extremely polyphagous species able to exploit several different food sources (e.g., Cryptocephalus fulvus feeds on several plant species belonging to at least 14 different families). To compare the structure and diversity of the microbiota of Chrysomelidae with different breadth of the trophic spectrum, the selected species were divided into trophic classes (generalist and specialist) depending on the number of host plant families: (i) specialist includes species feeding on a single plant family; (ii) generalist includes species exploiting more plant families.

Table 1

Information on the analysed samples.

Species	Diet	Subfamily	Date	Country	Latitude	Longitude
Altica oleracea	Generalist	Alticinae	08/02/2010	Italy	45.794 N	9.250 E
Chaetocnema hortensis	Generalist	Alticinae	08/02/2010	Italy	45.794 N	9.250 E
Crepidodera fulvicornis	Generalist	Alticinae	08/02/2010	Italy	45.794 N	9.250 E
Cassida inopinata	Specialist	Cassidinae	28/06/2009	Italy	44.517 N	8.819 E
Hypocassida subferruginea	Generalist	Cassidinae	14/07/2008	Italy	42.796 N	11.242 E
Dicladispa testacea	Specialist	Cassidinae	03/06/2011	Italy	44.195 N	8.281 E
Hispa atra	Generalist	Cassidinae	26/06/2011	France	42.512 N	2.124 E
Clytra quadripunctata	Generalist	Chryptocephalinae	17/07/2010	Italy	45.941 N	9.416 E
Cryptocephalus fulvus	Generalist	Chryptocephalinae	06/04/2010	Italy	40.855 N	12.956 E
Cryptocephalus loreyi	Generalist	Chryptocephalinae	17/05/2009	Italy	45.857 N	9.253 E
Cryptocephalus transcaucasicus	Generalist	Chryptocephalinae	25/07/2009	Italy	44.702 N	7.142 E
Labidostomis longimana	Generalist	Chryptocephalinae	13/07/2010	Italy	45.824 N	9.279 E
Pachybrachis exclusus	Specialist	Chryptocephalinae	21/06/2008	Italy	44.056 N	9.832 E
Smaragdina affinis	Generalist	Chryptocephalinae	17/05/2009	Italy	45.857 N	9.253 E
Calligrapha sp.	Unknown	Chrysomelinae	01/04/2017	USA	32.268 N	110.808 W
Chrysolina fastuosa	Generalist	Chrysomelinae	17/05/2009	Italy	45.857 N	9.253 E
Chrysomela saliceti	Specialist	Chrysomelinae	21/06/2011	Italy	44.456 N	9.823 E
Prasocuris phellandrii	Generalist	Chrysomelinae	24/04/2010	Italy	45.795 N	9.216 E
Timarcha tenebricosa	Generalist	Chrysomelinae	04/06/2012	France	43.847 N	6.518 E
Crioceris paracenthesis	Specialist	Criocerinae	04/06/2012	Italy	45.894 N	13.551 E
Lilioceris merdigera	Generalist	Criocerinae	21/05/2010	Italy	45.794 N	9.250 E
Donacia obscura	Specialist	Donacinae	11/05/2011	Italy	44.625 N	9.542 E
Plateumaris consimilis	Generalist	Donacinae	07/04/2017	Italy	45.794 N	9.250 E
Chrysocus asclepiadeus	Specialist	Eumolpinae	28/06/2010	Italy	45.831 N	9.286 E
Macrocoma henoni	Unknown	Eumolpinae	23/05/2013	Morocco	31.150 N	5.393 W
Exosoma thoracicum	Generalist	Galerucinae	11/06/2011	Turkey	37.232 N	27.611 E
Luperus longicornis	Generalist	Galerucinae	28/06/2010	Italy	45.831 N	9.286 E
Orsodacne cerasi	Generalist	Orsodacninae	28/05/2009	Italy	45.893 N	9.281 E
Orsodacne humeralis	Generalist	Orsodacninae	21/05/2011	Turkey	39.761 N	27.571 E
Zeugophora flavicollis	Specialist	Zeugophorinae	08/07/2011	Italy	45.943 N	9.410 E

Taxonomic (species, families), ecological (diet spectrum width, i.e., specialist or generalist) and collection (date, country, latitude, longitude) information are reported.

DNA extraction

DNA was extracted from the whole insect body using the classical phenol–chloroform methods (Doyle and Doyle, 1990) with the following modifications. First, 500 μl of 2% CTAB (2% CTAB, 0.2% ascorbic acid, 1.5% PVP, 1.4 mmol L−1 NaCl, 20 mmol L−1 EDTA and 100 mmol L−1 Tris–HCl, pH 8.0) was added to each sample. Tissues were then disrupted using glass beads (ø 0.1 mm) with the Precellys®24 homogenizer (Bertin Technologies, Montigny‐le‐Bretonneux, France) and incubated at 65°C for 15 min to inactivate nucleases. After centrifugation, the supernatant was incubated overnight with 20 μl of proteinase K (20 mg ml−1) at 56°C. To purify the DNA, two phenol–chloroform washes (phenol/chloroform/isoamyl alcohol, 25:24:1, pH 8.0) were performed. DNA was, then, precipitated after addition of 500 μl of isopropanol and incubation for 1 h. Pellet was washed twice with 70% ethanol and eluted in 40 μl of Ultrapure Water (Sigma‐Aldrich, Saint Louis, Missouri, USA). Qubit 4.0 fluorometer (Thermo Fisher Scientific) was used to determine the DNA concentrations. A DNA extraction blank, using the same extraction protocol and molecular biology grade water, was performed as control to monitor environmental contamination.

Library preparation and sequencing

Two regions of 16S rRNA gene (V1–V2 and V4) were sequenced by means of Ion Torrent platform (Life Technologies). PCR primers 27FYM (Frank et al., 2008) and 338R (Amann et al., 1990) were used to amplify V1–V2 region, while primers 515F (Caporaso et al., 2011), 802R (Claesson et al., 2009) and 806R (Caporaso et al., 2011) were used for V4 region, in two separated reactions. PCR primers were tailed with two different GC rich sequences enabling barcoding in a second amplification. The first PCR amplification of the V4 region has been performed as reported in Chouaia et al. (2019). The first V1–V2 PCR was performed in the same conditions as V4 following 34 cycles of 94°C for 15 s, 60°C for 15 s, 72°C for 15 s and a final extension of 72°C for 2 min. The second PCR amplification was performed in 25 μl volume containing 10 μl HotMasterMix 5Prime 2.5× (Quanta Bio), 1.25 μl EvaGreen™ 20× (Biotium), 1.5 μl barcoded primer (10 μM), 1 μl of the first PCR amplification with the following conditions: 8 cycles of 94°C for 10 s, 60°C for 10 s, 65°C for 40 s and a final extension of 72°C for 3 min. To control for bacterial contaminations, PCR amplifications of the V1–V2 and V4 regions were performed, as previously reported, using as template the DNA extraction blank (i.e., reagents of the used DNA extraction kit) and the PCR reagents. No amplicons were obtained by visualization on 1.5% agarose gel electrophoresis. Furthermore, real‐time PCRs were performed on the DNA extraction blank in order to monitor for contamination and select the appropriate number of the first PCR cycles to avoid the raise of the negative control curves. All the amplicons were checked for their quality and size by agarose gel electrophoresis, quantified with the Qubit™ dsDNA BR Assay Kit in the Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and pooled together in equimolar amounts. The library was purified running it in a precasted E‐Gel® SizeSelect™ (Invitrogen) agarose gel 2% and finally quality checked and quantified with High Sensitivity DNA reagents in the Agilent 2100 Bioanalyzer (Agilent Technologies). For sequencing the library was first subjected to emulsion PCR on the Ion OneTouch™ 2 system using the Ion PGM™ Template Hi‐Q OT2 View (Life Technologies) according to the manufacturer's instructions. Then ion sphere particles (ISPs) were enriched using the E/S module. Resultant live ISPs were loaded and sequenced on an Ion 316 chip (Life Technologies) in the Ion Torrent PGM System.

Bioinformatic analyses

The bioinformatic analyses were performed using the QIIME2 platform (Bolyen et al., 2019). The obtained raw reads for the two 16S rRNA gene regions (V1–V2 and V4) were denoised and taxonomically annotated separately. The DADA2 algorithm (Callahan et al., 2016) was used for denoising to obtain an estimation of the actual ASVs present (ASVs) using default parameters (e.g., trunQ = 2, maxEE = 2). The obtained ASVs have been taxonomically annotated with the fit‐classifier‐sklearn method (Pedregosa et al., 2011; Bokulich et al., 2018) using the release 138 of the SILVA database (Quast et al., 2012) as reference for sequences and taxonomy. The naïve Bayes classifiers were trained on the reference sequences trimmed to correspond to the amplified region. To obtain a common phylogeny for the ASVs from the two 16S rRNA regions, the SEPP technique (SATé‐enabled phylogenetic placement; Janssen et al., 2018) was applied to place the ASVs on a reference phylogeny based on the release 138 of the SILVA database (Quast et al., 2012). In specific cases (possible Enterobacteriaceae primary symbionts) the taxonomic annotation of the ASVs have been checked using the BLAST algorithm (Altschul et al., 1990) on the NCBI nt database and confirmed using phylogenetic tree inference.

To infer maximum‐likelihood trees for the phylogenetic placement of putative primary symbionts, 16S rRNA reference sequences for selected genera representative of the Enterobacteriaceae were downloaded from the NCBI nt database. Sequences were aligned using the mafft algorithm v.7.471 (Katoh and Standley, 2013) considering information on the secondary structures of the rRNA. The ML trees have been inferred with iq‐tree v.2.0.3 (Minh et al., 2020) using ModelFinder (Kalyaanamoorthy et al., 2017) to select the substitution model according with AIC (Akaike, 1973). Ten trees for each amplified region have been inferred to check for concordance of different runs. The same phylogenetic tree inference pipeline has been applied to manually aligned COI sequences of the 30 selected Chrysomelidae species (Magoga et al., 2018) (Supplementary Table 3) with the topology constrained to the one published in Nie et al. (2020).

The microbiota diversity analyses were performed with a sample size and coverage‐based integrations of interpolation (rarefaction) and extrapolation (prediction) of the Hill numbers (Hill, 1973; Alberdi and Gilbert, 2019; Roswell et al., 2021) using the R packages iNEXT and iNextPD (Chao et al., 2014, 2015; Hsieh et al., 2016). The computation of Hill numbers was performed for three increasing values of the order parameter q, corresponding to increasing weight on the species abundance (or any other taxonomic level considered) and also to different well‐known diversity and phylogenetic diversity indices: q = 0, counting mainly the rare species (those with low abundances), corresponds to richness (McIntosh, 1967) and Faith's Phylogenetic Diversity (Faith, 1992); q = 1, counting mainly the common species (those with medium‐high abundances), corresponds to the exponential of Shannon index (Shannon, 1948) and Allen's Phylogenetic entropy (Allen et al., 2009); q = 2, counting mainly the dominant species (those with very high abundances), corresponds to the inverse of Simpson index (Simpson, 1949) and Rao's quadratic entropy (Rao, 1982). This explicit parametrization is particularly useful to test our hypothesis, since we can assume that symbionts with a functional role are present at high abundances (q = 1, q = 2) while the bacteria acquired from the environment have usually low abundances, so can be considered rare species (q = 0). The hilldiv R package (Alberdi and Gilbert, 2021) was used to partition the diversity in its components: α‐diversity (diversity at the sample level), γ‐diversity (total diversity in the selected group) and β‐diversity (γ‐diversity/α‐diversity, corresponding to the among sample component of the total diversity). Comparing diversity partitioning between the two trophic groups considered in this study (specialists and generalists) allows us to understand which component is most influential in determining the different diversity estimates. The R package phytools (Revell, 2012) was used for ML ancestral state reconstruction (Revell, 2013) of the microbiota diversity estimates along the insect phylogenetic tree and to compute phylogenetic signals (Blomberg et al., 2003; Ives et al., 2007).

Supplementary Fig. 1. Bacteria abundance in single marker datasets. Heatmap representing the abundance of bacterial taxa (classes, families, genera) present in the single marker datasets (V4 and V1‐V2). In the genera heatmap only the 50 most abundant genera are shown. Colour intensity is proportional to the normalized relative abundance of the bacterial taxa.

Click here for additional data file.

Supplementary Fig. 2. Microbiota diversity estimates inferred on the total dataset (V1‐V2 and V4), V1‐V2 and V4 regions of the 16S rRNA. Sample‐based rarefaction/extrapolation curves of the Hill numbers estimated for three values of the order parameter (q = 0, q = 1, q = 2). The x‐axis represents increasing sampling and the y‐axis represents the Hill number estimates, 95% confidence interval is also reported. As reported in the legend, colours correspond to the trophic category (specialist or generalist) and line type to the methodological approach (interpolation or extrapolation). a) Global dataset (V1‐V2 and V4 regions of the16S rRNA). b) V1‐V2 region of the 16S rRNA. c) V4 region of the 16S rRNA.

Click here for additional data file.

Supplementary Fig. 3. Maximum likelihood phylogenetic trees. Sequences obtained from the NCBI database report the accession numbers while sequences produced in this study are highlighted in bold. a) Tree obtained from sequences of the V1‐V2 region of the 16S rRNA. b) Tree obtained from sequences of the V4 region of the 16S rRNA.

Click here for additional data file.

Supplementary Fig. 4. Microbiota diversity of specialist and generalist Chrysomelidae defined using the plant taxonomic level of genus (specialists feed on plants all belonging to the same genus, generalists feed on plants belonging to different genera). Coverage based rarefaction/extrapolation curves of the Hill numbers estimated for three values of the order parameter (q = 0, q = 1, q = 2). The x‐axis represents the coverage (that estimates the completeness of the sampling) and the y‐axis represents the Hill number estimates, 95% confidence interval is also reported. As reported in the legend, colours correspond to the trophic category (specialist or generalist) and line type to the methodological approach (interpolation or extrapolation).

Click here for additional data file.

Supplementary Table 1. Primary symbionts relative abundance.

Supplementary Table 2. Blast search results on NCBI nt database.

Supplementary Table 3. Accession numbers of COI sequences.

Click here for additional data file.