Samples

We sampled several species of leaf beetles specialized on alder, specifically found on speckled alder, A. incana ssp. rugosa, in three localities in the area between the counties of Franklin and Essex, in the Whiteface Mountain region, upstate New York (Fig. 1; Table 1). In particular, we obtained samples of Calligrapha alni Schaeffer, C. confluens Schaeffer, Chrysomela mainensis ssp. mainensis Bechyné, and Altica ambiens ssp. alni LeConte. None of these species has been tested before for the presence of endosymbionts in any part of their range. Moreover, in two of the localities, we also obtained samples of leaves from the trees where the beetles were captured to investigate the diversity of bacteria on their surface. In the case of the insects, we avoided manipulation and they were transferred with forceps from the leaf where they were sitting to individual tubes containing absolute ethanol; in turn, fragments from central parts of some leaves were cut with scissors, while holding the fragment with clean forceps, and also directly transferred into a vial with ethanol. The sample was completed for exploratory purposes with one specimen of Calligrapha philadelphica (Linnaeus) from one of the localities (Table 1), a sympatric species but with a completely different trophic selection, since it is exclusively associated with dogwood, Cornus spp. (Cornaceae).

Fig. 1.

Geographic location of the study sites. Locality codes and collection data are the same given inTable 1.

16S rRNA Gene Metabarcoding Library Preparation and Sequencing

Owing to the effects of ethanol, including hardening and loss of structure and color of internal organs, it was not possible to recognize features of the internal anatomy of the leaf beetles. Thus, we dissected under the microscope and in sterile conditions tissue fragments of the abdomen of each specimen, hoping to extract enough tissue typically carrying endosymbiont colonies, preferentially the gonads (Pietri et al. 2016). The dissection was facilitated soaking the tissues in PBS 1× or the first buffer (AHL) of the DNA extraction kit. In particular, microbial genomic DNA was obtained from each leaf beetle sample individually using the QIAamp DNA Microbiome Kit (Qiagen Iberia S.L., Barcelona, Spain) and following the manufacturer's instructions. In the case of the leaf samples, we used approximately 1 cm² of preserved leaf fragments from Hough Brook and the Flume Trail system, which were subject to the same bacterial DNA extraction protocol as the beetle tissues.

The 16S rRNA gene metabarcoding libraries targeted the V3 and V4 regions of the large ribosomal RNA subunit and using the primers of Klindworth et al. (2013) with Illumina adapter overhang nucleotide sequences (in italics) as follows: Forward Primer, 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG; and Reverse Primer, 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC. The PCR used some 10–15 ng of microbial DNA (except in the negative control), 0.2 µM final concentration of each primer, and the KAPA HiFi HotStart ReadyMix (Roche Diagnostics, Rotkreuz, Switzerland), with a standard PCR cycling protocol: 95°C for 3 min, followed by 25 cycles of 95°C for 30 s, 55°C for 30 s, 72°C for 30 s, and a final extension step at 72°C for 5 min. PCR products were purified and cleaned off primers and primer dimers with AMPure XP beads (Beckman Coulter, Indianapolis, IN), and 5 µl of cleaned PCR products used as template for an indexing PCR with the Nextera XT Index kit (Illumina, San Diego, CA) to attach dual indices and Illumina sequencing adapters to each library. Indexed libraries were cleaned again with the AMPure XP beads and 1 µl of each library checked for quality and size on a Bioanalyzer DNA 1000 chip (Agilent Technologies, Santa Clara, CA). Libraries were normalized to a final concentration of 4 nM with 10 mM Tris pH 8.5, and 5 µl of each indexed library pooled for a single MiSeq run. Prior to sequencing, DNA of pooled libraries was denatured with 0.2 N NaOH, the amplicon library diluted with HT1 hybridization buffer according to Illumina protocols and combined with 5% of the PhiX Control v3 library. The pooled sample of 19 libraries was loaded on a MiSeq for a 2 × 300 bp paired-end sequencing run in the Genomics Core facility at the University Pompeu Fabra (Barcelona, Spain). Haplotypes discussed in the text have been deposited in the European Nucleotide Archive (EMBL-EBI, Hinxton) under sequence accession numbers LR877032–LR877144.

Table 1.

Sampling localities, species studied, and number of individuals

Microbiome Characterization

The sequencing procedure generated paired FASTQ files for each library that were sorted based on their indices, excluding all individual sequences with quality and indexing issues. Amplicon 16S rRNA gene primers and adapters were removed from each sequence after this initial filtering using the Python tool Cutadapt 2.6 (Martin 2011). Trimming ends of low-quality, fragment length filtering (truncLen=c(225,205), maxEE=c(2,2), truncQ=2), merging of paired ends, removal of length outliers by retaining fragments of 400–414 nucleotides, and removal of chimeras (based on the consensus method) were done using the R package DADA2 v1.14.0, which uses a procedure that models and corrects sequenced amplicon errors in Illumina (Callahan et al. 2016). The taxonomic identification of the resulting data sets was done based on the March 2018 update of the Silva 132 reference database, the one currently with highest taxonomic coverage (Quast et al. 2013), and the IDTAXA algorithm for classification (Murali et al. 2018), as implemented in the R package DECIPHER v2.14.0 (Wright 2016). Among the advantages of DECIPHER, it handles sequence data as a database, instead of reading sequences to the computer memory, making it possible to handle high-throughput sequence data efficiently even with computers that are not high-performance (see below). In turn, IDTAXA exploits machine-learning algorithms which reduce misclassification errors by identifying sequence diagnostic features of specific ranks in the training set, which are also statistically rare in others. These traits are used to reassess the success in the classification of the actual training set and assign confidence values to the classification procedure by traversing the nodes of the classification hierarchy from the root to the tips. The same procedure is applied to query sequences by requiring strict (bootstrap support = 98%) associations of a query within a particular level of the classification to yield that result as the final taxonomic assignment.

Phylogenetic Identification Pipeline

The previous taxonomic identification pipeline allowed sorting the sequences in major taxonomic groups and singling out the generally abundant endosymbiont sequences. However, in order to find finer taxonomic structure in the endosymbiont data, we used an adaptation of our automated DNA-based phylogenetic identification pipeline, BAGpipe (Papadopoulou et al. 2015), on the full availability of bacterial sequences from GenBank. This pipeline was originally designed to study psbA-trnH cpDNA plant sequences, but it can be relatively easily adapted to other markers and other organisms (Papadopoulou et al. 2015). The first stage of the pipeline consists of creating a reference database from GenBank data for the marker of interest, matching orientation and length of the studied barcode, and removing taxonomic redundancy, while retaining relevant taxonomic information for downstream identification and filtering steps. We followed the pipeline instructions for database construction with some modifications and all the analyses were run on an iMac OS X 10.7.5, with 2.7 GHz Intel Core i5 processor and two 2 GB memory modules (1333 MHz DDR3). The 387 flatfiles of the latest (November 2019) bacteria GenBank release (gbbct) and the associated GenBank taxonomy (NCBI taxon number for bacteria = 2) were downloaded to automatically parse all the available V3–V4 16S rRNA gene fragments. Identification of the homologous fragment of interest was based on an anchor query text file with a selection of V3–V4 16S rRNA gene sequences from our own libraries, as sequenced and identified in the previous steps. Some steps for the preparation of the reference database, namely reorientation and trimming of sequences using USEARCH (steps 1.6 to 1.8 of the original pipeline description), required splitting the original files that included all the available sequences into smaller files due to processing memory limitations in our setup, possibly because of many entries corresponding to complete genomes, thus also including the V3–V4 16S rRNA gene region of interest. In particular, in this implementation of the pipeline, we used USEARCH v5.2.32 (Edgar 2010) because the version recommended in the description of the pipeline (v4.66) was deprecated and no longer available.

The second stage of the pipeline consists of the identification of query sequences by comparison with the previous reference database based on either phenetic, phylogenetic, or both criteria, by reducing the problem to manageable groups of similar sequences in the query group. In our case, we split the endosymbiont sequences and homologous sequences from the reference database in 90% similarity clusters, which fundamentally distinguished between two types of endosymbiont, and submitted them to phylogenetic clustering and identification. Query and reference sequences were aligned with the G-INS-i algorithm in MAFFT 7 (Katoh and Standley 2013), and the resulting matrices used for maximum likelihood phylogenetic inference using RAxML 7.2.8 (Stamatakis 2006) implementing a GTR+CAT model on a heuristic hill-climbing search based on 100 random most parsimonious initial trees, and assessing node support based on 100 bootstrap pseudoreplicates. Also, the relationships among the empirically obtained haplotypes of the two endosymbionts present in our sample were investigated using statistical parsimony networks based on the software TCS 1.21 (Clement et al. 2000).

Characteristics and Quality of Library Reads

16S rRNA gene amplification and library preparation using the multi-tagged primers yielded in every case fragments of approximately the expected size (592–607 bp) and relatively low, but usable concentrations (2.74–16.31 ng/µl) (Table 2). Index-based sorting produced on average nearly 875,000 reads per amplicon library, ranging between a low ∼565,000 reads in one sample of Ca. alni, and the expected ∼1,110,000 reads in two samples of Ch. mainensis (Table 2). Filtering based on sequence end-quality threshold, denoising based on the explicit error rate models of the procedure that we followed, merging of paired ends, and removing length outliers retained 43.9–68.0% of the initial reads. Finally, 0.33% of the surviving reads were automatically recognized as chimeras and also excluded from subsequent analyses (Table 2).

Table 2.

Parameters of 16S rDNA library synthesis based on the Bioanalyzer results and of the sequence-filtering procedure with DADA2 (Callahan et al. 2016)

Composition of Libraries

One striking result after the initial taxonomic assignment of all sequences based on the reference database Silva 132 was finding that more than half of the libraries were overwhelmed with contaminant plant mitochondrial, but mainly plastid V3–V4 DNA sequences. Interestingly, this problem did not affect samples randomly, and it fundamentally showed in all Calligrapha and leaf samples, while most samples of the other two leaf beetle species included normal values between 0.2 and 66.7% of plastid and mitochondrial sequences (Table 2).

Up to 36 sequences (7.47% of all bacterial reads) were not assigned any taxonomic rank or yielded very high ranks in the taxonomic hierarchy (e.g., Gammaproteobacteria). Most (98.45%) of these unidentified reads actually corresponded to a unique sequence obtained from one of the libraries of A. ambiens (and just eight reads from one of the alder leaf samples), which showed maximum similarity (approximately 94%) with environmental sequences of saline habitats inhabited by natural populations of Artemia in Israel (Tkavc et al. 2011). Of the remaining 35 unidentified sequences, 15 (197 reads) were recognized in the alignments or with additional BLAST searches as chimeras of highly divergent taxa, and eight (89 reads) were confirmed as plastid sequences affected by a particular sequencing problem: stretches of 50–100 bp where all G residues were replaced by T (or C to A, in the reverse strand). More worrying, although always relatively considering their minimal representation in the sample (61 reads), were five sequences which were attributed family ranks in the Lactobacillaceae (Bacilli) and the Porticoccaceae (Gammaproteobacteria). Multiple sequence alignments of all ingroup 16S rRNA gene sequences and BLAST searches helped recognizing them as chimeras between sequences of Wolbachia and of bacteria of the groups mentioned.

After excluding plant contaminant sequences and the obvious artifacts mentioned above, we found evidence for 115 different microbial 16S rRNA gene sequences in the whole data set (42 operational taxonomic units [OTUs] at a 97% identity threshold), of which 58 (corresponding to 91.51% of all bacterial reads) were unambiguously identified as Rickettsiales of the genera Wolbachia and Rickettsia. Apart from the abundant Wolbachia and Rickettsia, the highest diversity of bacterial identifications corresponded to the samples of Ch. mainensis. In this species, Actinobacteria of the Frankiales (Nakamurella), Micrococcales (Brachybacterium and Kocuria), and Propionibacteriales (Cutibacterium), and Alphaproteobacteria of the Acetobacterales (e.g., Acetobacter) and of the Rhizobiales (Afipia and Rhodopseudomonas) dominated the sample (Fig. 2;  Supp Table S1 [online only]). The next most diverse sample, albeit deduced from a reduced number of observations and again by excluding Wolbachia, corresponded to the bacteria characterized from Calligrapha confluens. In this case, we also found evidence for associations with Kocuria and Cutibacterium among the Actinobacteria, Rhizobiales of three different families (A0839, Hyphomicrobiaceae and Rhizobiaceae), as well as Clostridia and different types of Gammaproteobacteria (Fig. 2). In the case of Altica, and excluding the Rickettsiales, the sample was dominated by Actinobacteria of the Mycobacteriaceae (Mycobacterium), Micrococcaceae, and Propionibacteriaceae (Cutibacterium), but it also had evidences for Cyanobacteria and Proteobacteria of the Rhizobiales (Methylobacterium). Finally, Alnus leaf surfaces yielded a rather diverse set of inferences, mainly Alphaproteobacteria of the Acetobacterales (Acetobacteraceae, including Acidiphilium and Endobacter) and of the Rhizobiales (Beijerinckiaceae, including Methylobacterium and the unidentified genus 1174-901-12).

Diversity of Wolbachia and Rickettsia

The sample yielded 36 sequence variants unambiguously assigned to Wolbachia and 22 to Rickettsia, with the former retrieved from all beetle species, and the latter predominantly in A. ambiens, but with some reads in the libraries from Ch. mainensis and Ca. alni. Interestingly, for one of the leaf samples, most of the bacteria characterized also belonged to Wolbachia, and specifically with the same most abundant 16S rRNA haplotype as retrieved from the beetle samples. Despite the relatively large number of 16S rRNA haplotypes, four haplotypes (W1–W4) represented 98.5% of all Wolbachia sequences, and a single haplotype 95.9% of the Rickettsia contigs. Most parsimonious networks relating these haplotypes revealed a number of interesting patterns (Fig. 3). In the case of Rickettsia (Fig. 3a), the dominant haplotype (R1) was the only one also found, if at very low frequencies, in Ch. mainensis and Ca. alni, and this haplotype was the only one shared between the two specimens of A. ambiens tested, which otherwise had their discrete sets of haplotypes (R2–R4 in one individual, and all the others in the other). For Wolbachia, 16S rRNA genetic diversity was higher than for Rickettsia and haplotypes belonged to two main groups separated by 10 mutational steps, one with haplotypes W4 (one of the most abundant in the sample) and W36, exclusively found in A. ambiens, and a group with all the other haplotypes found in all species tested (Fig. 3b). The last group showed a central core with the other three dominant haplotypes from which radiated all others in star-shaped patterns, particularly from the most abundant W1. Each of these haplotypes defined nonetheless three domains, one derived from W3 and nearly exclusive of A. ambiens, one derived from W1 for Ch. mainensis, and one derived from W2, dominated by Ch. mainensis, but with haplotypes W28 and W29 exclusively found in Calligrapha spp. Interestingly and for the reasons that will be elaborated in the discussion, the latter haplotypes were the only ones recovered from the library of Ca. philadelphica.

Fig. 2.

Taxonomic diversity and abundance of nonendosymbiotic bacteria deduced from 16S rRNA haplotypes from the libraries of four species of the study. The samples of Calligrapha alni and Ca. philadelphica only yielded sequences of Rickettsia and Wolbachia and are not shown here. Frequency axes apply to Chrysomela mainensis and Altica ambiens (left) and to Ca. confluens and Alnus incana (right).

The 16S rRNA gene phylogeny of Rickettsia showed very little resolution, but four main groups were distinguishable: one of them only characterized in the whitefly B. tabaci and other hemipterans; one mainly found in weevils; one in diverse organisms, but predominantly water beetles and leeches; and a last one, highly diverse, poorly structured and characterized from all sorts of insects and arachnids ( Supp Fig. S1 [online only]). The Rickettsia haplotypes characterized from A. ambiens and in lower frequencies in two other species of leaf beetles belong in the latter group. The only Rickettsia haplotype in our sample that was already present in the database was the most abundant R1, identical to sequences previously characterized in mites, lacewings, different genera of true bugs, fleas, and all major groups of holometabolous insects, including several groups of beetles. This same 16S rRNA haplotype was the one characterized from the type strains of R. bellii and R. australis.

Fig. 3.

Unrooted statistic parsimony networks of 16S rRNA haplotypes of Rickettsia (a) and Wolbachia (b) obtained for four species of alder specialists of the Chrysomelidae. Each colored node represents an observed haplotype and its code and absolute frequency are given, and the size of the node is proportional to the order of magnitude of this frequency. Edges represent mutational steps and small black circles missing intermediate states.

The diversity and phylogenetic structure of the 16S rRNA gene for Wolbachia as deduced from public sequence databases was higher than for Rickettsia, but still not enough for high resolution of supported phylogenetic groups ( Supp Fig. S2 [online only]). But relevant for the purposes of the current study, this taxonomic marker allowed discriminating between two supergroup clusters of the bacterium, including an assemblage of representatives in supergroups A, B, and E, where the haplotypes W4 and W36 exclusive of A. ambiens were found, and another assemblage with representatives of supergroups C, D, and F where all the other Wolbachia haplotypes found in this study clustered. Haplotypes W4 and W36 appeared in a large, poorly resolved group of haplotypes described from numerous insect, arachnid, and isopod groups, and specifically recognized as members of supergroup B in Wolbachia (e.g., Mateos et al. 2006). The remaining haplotypes appeared in a lineage within the supergroup F of Wolbachia with many endosymbionts predominantly characterized from ants and beetles, including weevils and other leaf beetles in the genera Basilepta, Galeruca, and Lema.

Technical Issues of Metabarcoding and the Reliability of Inferences

Before extracting the main conclusions from our study, we consider that it will be useful for other researches working in similar systems to be aware of and have some context for a number of artifacts that plagued our results. The problems that we detected in our study belong to two main, well-known drawbacks of metabarcoding studies (Kunin et al. 2010, Pinto and Raskin 2012, Tamošiūnė et al. 2019), including 1) artifacts derived from PCR, and 2) sequencing problems, with some subsequent issues in automatic identification procedures. A third type of problem, in part related to the PCR-based approach as well, was the amplification of nontarget sequences, which will be discussed in the next section.

Among the common problems of PCR-based approaches to metabarcoding, we could not assess here biases in species detection due to primer mismatches or bacterial contamination (e.g., von Wintzingerode et al. 1997, Hong et al. 2009). However, in the specific case of the study of Wolbachia and Rickettsia sequences, their analysis using genealogical approaches revealed patterns consistent with at least a certain proportion of PCR errors being resurfaced by the power of the sequencing technology used here. Theory predicts that PCR produces erroneous sequence variants due to point mutations, with their frequencies diminishing as they accumulate more changes (Weiss and von Haeseler 1995). These PCR errors typically result in star-like tree topologies with the target sequence as a more abundantly represented core and de novo sequences resulting from polymerase mistakes radiating from it, most of them by a single mutation (Smart et al. 2019). This is exactly the pattern that we retrieved, suggesting that a great deal of lower frequency haplotypes may indeed represent PCR artifacts and not true diversity of Wolbachia and Rickettsia in the sample. This idea is strengthened by the fact that these putative artificial sequences are absent in public sequence databases (with the exception of W28, common in weevils and also found in one leaf beetle; Toju et al. 2013, Jeong et al. 2019). Conversely, our core haplotypes have been characterized already in other organisms and are shared among species in our sample (W29 is also shared between three different Calligrapha species, thus suggesting that it may be real too). Focusing on the most abundant core haplotypes and not on the whole diversity revealed by the ultrasequencing approach used here, even though some of it may represent true endosymbiont diversity, makes our conclusions conservative and equally valid.

Some of the point mutations discussed above may be also the result of wrong base calling in the sequencing stage, rather than PCR artifacts. However, the most obvious sequencing problem that we detected, if for a reduced number of sequences, was the complete transformation of G residues in the original sequence to T residues in the deduced sequence for a relatively long stretch of nucleotides (and mirrored by changes from Cs to As at the opposite end of the sequence, thus in the paired read). The effect was observed in one Wolbachia, one Rickettsia, and several plastid sequences. Given that this specific problem affected both paired ends of the same sequence, we interpret it as a topological artifact, a local saturation with T fluorophore that outshined the G fluorophore, as their emission spectra overlap in the 4-channel technology used by MiSeq systems (fluorophore cross talk; Sheikh and Erlich 2012). This effect lasted until the excess of T was successfully washed away, and sequencing resumed normally. What is important, however, is that these sequences, even though they had a very low representation (0.005% of reads), passed the pipeline quality filters and reached the identification steps of the analysis, contributing 10 haplotypes (7%) and a number of unidentified OTUs to the global results.

Another problem derived from PCR was the formation of chimeras (Haas et al. 2011), which can be approached analytically in most cases, although it was revealed here because of some limitations in the quality control of library composition. We had evidence that a number of objective chimeras survived the quality-filtering steps of the identification pipeline, and only a meticulous analysis of individual 16S rRNA sequences in multiple sequence alignments, phylogenetic trees, and by BLAST searches against public sequence databases allowed us to effectively remove the most obvious cases. Using this approach we could filter out a number of obvious chimeric sequences (0.008% of reads; 10.5% of haplotypes) and identify their potential sources. Interestingly, even though most of these artifacts were not identified below the rank of order, in some cases, when most of the chimera corresponded to a particular identifiable taxon, they yielded a positive taxonomic identification at least at family level. As with the T-saturation problem mentioned above, only a manual examination of data made it possible to detect these sequences, often with partial Wolbachia sequences, and remove them from subsequent interpretations of data.

The Challenge of Metabarcoding Libraries of Herbivorous Insects

Although not a technical problem per se, the main setback that we found in our study was the hyperabundance of contaminant plastid sequences, which nearly obliterated target bacterial sequences in some libraries. Finding these specific plant contaminants is common and a permanent source of concern in metabarcoding studies of microorganisms from certain environments (e.g., Sakai and Ikenaga 2013). For example, metabarcoding studies of plant microbiomes, research of soil samples and of herbivorous insects have shown that organellar sequences can account for more than 80% of data retrieved from these samples (as reported by Lundberg et al. 2013). Conversely, studies on nonphytophagous beetles show that the proportion of these contaminants stays very low (Hammer et al. 2015). In our case, working with abdominal contents of herbivorous beetles and plant tissue, finding them was certainly expected. However, with a 99.9% rate of retrieval of these contaminant sequences in a number of samples, this surpassed any original expectation. Similar studies focusing on herbivore beetles typically report that chloroplast (and/or mitochondrion) sequences were detected and removed. However, the figures to assess the importance of these findings are often not shown or are buried in supplements (Hanshew et al. 2013). The few explicit data that we could find show that these figures are erratic, even within a single study, and it is difficult to recognize a pattern. For example, Hammer et al. (2015) found around 20% of chloroplast sequences in samples of phytophagous Epilachna ladybugs; Kelley and Dobler (2011) reported about half the sequences of their Longitarsus flea beetle microbiome libraries being of chloroplastic origin; or Blankenchip et al. (2018), in a study of microbial profiles in Cephaloleia leafminers, found up to 93.6% nontarget sequences in one of their libraries, and 60.6% (SD = 19.2%) on average in their whole data set. Such lack of consistency possibly responds to idiosyncratic factors from each sample and/or study, in part related to practical choices (e.g., tissues used for library preparation), incomparability of approaches (e.g., PCR primers), metabolic state of the individual (e.g., food quantity and time since last meal), or ecological determinants (e.g., host plant species), among others.

Beckers et al. (2016) reviewed this problem and the strategies that researchers use to avoid it. These include 1) DNA extraction protocols that minimize organellar DNA co-extraction, 2) use of blocking agents to prevent PCR amplification of specific targets, or 3) the method of choice in most studies, based on mismatch primers specifically designed to avoid the amplification of organellar genes. We opted for the latter strategy, using V3–V4 mismatch primers supposedly excluding eukaryotic versions of the 16S rRNA gene, at least in environments with low plastid input (Klindworth et al. 2013, Beckers et al. 2016). Yet, this strategy proved inefficient in our hands, at least for some of the samples, specifically all Calligrapha individuals tested and one of Ch. mainensis. In those cases where we obtained a high representation of plastid genes, it is likely that they were preferentially amplified because they were already in a much higher proportion relative to bacterial genes in the extracted DNA (Sakai and Ikenaga 2013). In retrospective, it seems possible that our dissections of abdominal tissue in the larger Calligrapha beetles ended up taking gut tissue, and unfortunately also gut contents, e.g., recently consumed and partially digested plant tissue. A lesson directly derived from the study of Fitzpatrick et al. (2018) on the efficiency of peptide nucleic acid clamps to prevent nontarget coamplification depending on the organelle gene sequence is that the efficiency of mismatch primers may also depend on the plant species. We cannot discard that the unintended presence of plant tissue coupled with a reduced mismatch efficiency in alder contributed to exacerbate the nontarget amplification problem.

Regardless of the actual explanation for such a high representation of nontarget sequences in some of our libraries, our main concern was whether the idiosyncratic composition of these libraries could diminish the reliability of the comparatively few endosymbiont (or other bacterial) inferences retrieved. Several lines of evidence suggest that these sequences may indeed represent trustworthy associations. One of the Ch. mainensis libraries showed the same nontarget near exclusivity as seen in all samples of Calligrapha, but it also showed evidences for a relatively small number of Wolbachia sequences. These contigs belonged without exception to haplotype W1, the most abundant by far in all the other Chrysomela samples. In turn, Calligrapha libraries had a low representation of Wolbachia and/or Rickettsia sequences, but they were either identical to the most abundant sequences characterized in other species or were shared between Calligrapha species, enormously reducing the odds of being an artifact. The same would be true for the bacteria characterized on the leaf samples, particularly those of Wolbachia, again belonging to the W1 haplotype, the most abundant in the entire sample.

Fig. 4.

Schematic representation of endosymbiont diversity and shared associations of Rickettsia (R haplotypes) and Wolbachia (W haplotypes) among alder specialists of the study. Haplotype W1 of Wolbachia is shared among all species, including samples of the host plant, Alnus incana, and haplotype R1 of Rickettsia was detected in the three of the studied genera of Chrysomelidae. [Pictures used with permission from the authors: see Acknowledgments.]

Implications of Endosymbiont Diversity and Their Distribution

Some remarkable results from our study include 1) finding that all samples, except for two of the three leaf extracts, had evidence for Wolbachia, for Rickettsia, or for both, 2) that these bacteria were clearly dominant in all cases, and 3) that more than one type of endosymbiont were present in every sample (Fig. 4). The rate of infection with Wolbachia in Coleoptera varies considerably, depending on the host taxonomy, with Chrysomelidae showing intermediate values (Kajtoch and Kotásková 2017). However, our study precisely focused on two genera, Altica and Calligrapha, which have been shown to have particularly high infection rates with Wolbachia among the leaf beetle groups tested (Jäckel et al. 2013, Gómez-Zurita 2019). Rickettsia has received much less attention than Wolbachia in studies of associations of bacteria with Coleoptera. Yet, although it is less prevalent than Wolbachia, it is not rare to find this association too (Kolasa et al. 2019). The fact that these endosymbiont Rickettsiales are prevalent in the specimens studied could explain in part the low diversity observed for other bacteria (Kolasa et al. 2019). Also in support for this idea, Blankenchip et al. (2018) empirically found that beetles infected by Rickettsia and Rickettsiella had lower levels of other bacterial groups, while beetle species free from Rickettsia presented microbiomes more diverse. This particular effect of endosymbionts modulating the diversity of other bacterial associations could be coupled with several other characteristics of our system that have been related with lower numbers of host associations with bacteria. For example, gut bacterial diversity has been explained by a number of factors including the taxonomic group (Yun et al. 2014), and in this respect, leaf beetles seem to have typically lower diversity in their microbial profiles compared to other beetle groups (Kolasa et al. 2019). Also, generalist species have been proposed to harbor more diverse microbiota, perhaps related to their exposure to a wider range of bacterial sources (Engel and Moran 2013, Hansen and Moran 2014; but see Blankenchip et al. 2018). However, this latter explanation may be questionable, since individuals of generalist insect species are nonetheless most likely exposed to a single or very few potential host plants through their lives. In any case, the species studied here, known alder specialists (Clark et al. 2004), have necessarily restricted diets and exposure.

Lessons About Endosymbiont Transmission Among Alder Specialists

The presence of more than one endosymbiont in each sample and information on their relative proportions represent the findings with more relevant implications relative to the mode of transmission of these bacteria. It appears that each species has their own typical endosymbiont strains of either Wolbachia and/or Rickettsia: W1 and W2 Wolbachia strains in Chrysomela; W28 and W29 Wolbachia strains in Calligrapha; and W3 and W4 Wolbachia and R1 Rickettsia strains in Altica. Given the taxonomic conservatism and their prevalence in every case, we can assume that these associations represent the stable, vertically inherited germline infections for these individuals or populations. Differences in the life histories of these insects and the ways they relate with the environment could explain the patterns of endosymbiont exclusivity and also potential sources of transmission from external sources in the trophic community. In principle, one could speculate that these differences would have a major impact if they affected larvae, when individuals are first exposed to external infections. However, the most notable difference among these species has to do with pupation, occurring underground in the case of Altica and Calligrapha, and on the host leaves in the case of Chrysomela. Otherwise, these beetles lay eggs on the host plant; develop as larvae feeding on leaves (Altica are rather an exception among flea beetles, which usually develop as larvae in the soil; Doguet 1994); and live as adults associated with the host plant (Woods 1917, Whitehead 1919, Brown 1956). There is not much information on endosymbiont turnover across these developmental stages in insects, but it is possible that contrary to gut bacteria (Yun et al. 2014), pupation and the drastic disintegration of larval tissues do not necessarily involve a break in the association with endosymbionts, which are carried along from larval to adult stages (e.g., Goto et al. 2006, Zhukova et al. 2017). If this were the case for these leaf beetles, life-history similarities would rule out most species-specific external sources of endosymbiont transmission (except, perhaps, the contribution of parasitoids; Ahmed et al. 2015).

Endosymbiont exclusivity contradicts in part our initial expectation that a tight ecological network like the one portrayed in this study, where different species can be found contemporaneously on a same plant, may have higher chances to share their endosymbionts as well. These results mirror those obtained for other phytophagous beetles in equally close trophic communities, competing for identical resources (Merville et al. 2013). But, what do the lower frequencies of these Wolbachia and Rickettsia haplotypes represent in other species? Considering that these endosymbionts could reflect either stable infections or transient associations with the host (e.g., Chiel et al. 2009, Chrostek et al. 2017, Joubert and O'Neill 2017), can we discard their acquisition from the environment and challenge the hypothesis that common food choice favors sharing of endosymbionts? The strain represented by haplotype W1 in Wolbachia is the most informative to advance an answer to these questions thanks to a couple of circumstantial evidences: 1) this strain of Wolbachia was also confidently found on a leaf extract, while it was not supposed to be there because at the time of collection or during processing it had not been in touch with any of the leaf beetles; and 2) the only sample outside of the chrysomelid-alder system, i.e., the sample of the dogwood specialist Ca. philadelphica, was also the only one without W1 Wolbachia, although it objectively shared the W28/W29 strains, considered vertically transmitted and stable associations with the other two species of Calligrapha in the study. The latter is interesting, since it relates this finding with the other only study on Calligrapha endosymbionts, because we know that individuals of Ca. philadelphica in the locality of the study and in the general area are infected by two Wolbachia strains that were named wCalA1 and wCalC2, with wCalA1 shared with at least three other species, two of them specialized on yet a different host plant, willows (Gómez-Zurita 2019). These data are compatible with the hypothesis that host plants favor shared endosymbiont associations, suggesting a number of interesting ideas, worth testing with a more extensive study. These include, for example, that W1 in this system is in a stable association with Ch. mainensis based on its high representation and in all individuals analyzed. Also, that these bacteria are transitorily present in the environment (e.g., Alnus leaves), most likely propagated from their usual host. And, finally, that other herbivore beetle species living in the same environment get exposed to and introduce in their bodies these abundant bacteria, simply by feeding. Similar statements could be proposed for the R1 variant of Rickettsia, characteristic of A. ambiens and found at low frequencies in Ca. alni and Ch. mainensis; however, we prefer to be more cautious about interpretations on this endosymbiont, since we did not find it on Alnus leaves.

Based on metabarcoding studies, particularly using powerful next-generation sequencing (NGS) methodologies and a low-resolution taxonomic marker such as the 16S rRNA gene, we cannot know whether the secondary associations that we detected also respond to secondary infections of new hosts, thus accomplished horizontal transfers, or just a transient stage in the gut of the studied individuals after they fed upon contaminated food (e.g., Chiel et al. 2009). However, even if we could prove that endosymbionts like the one represented by the W1 strain of Wolbachia in our system were not transferred horizontally and infected alternative hosts, that it was a transient association, data from other studies support that this mechanism is possible (e.g., Mitsuhashi et al. 2002, Sintupachee et al. 2006, Fenton et al. 2011, Caspi-Fluger et al. 2012, DeLay et al. 2012, Gonella et al. 2015, Li et al. 2016, Kolasa et al. 2017). In this context, we can still contribute additional understanding about this transmission mechanism. At least the first necessary stages in the chain of events leading to horizontal transfer via food, i.e., the transmission from an original host to the environment and the acquisition from this external source, a host plant in this particular case, by a new insect host is certainly a possibility (Chrostek et al. 2017). Relative to the presence of the bacteria on the plant, we still favor a viewpoint that considers this to be accidental and temporary, e.g., through feces or laid eggs. However, new evidence suggests that endosymbionts may have some potential to survive in environments outside of the insect host, including plant tissues (Caspi-Fluger et al. 2012, Gonella et al. 2015, Li et al. 2016). In any case, this possibility only increases the chances to acquire endosymbionts through ingestion even from phylogenetically distant hosts that use the same host plant.

Conclusions

The accessibility, abundance, and taxonomic diversity of insects associated with Alnus and the high natural prevalence of Wolbachia and Rickettsia in species assemblages like the one studied here contribute to this trophic network being a suitable system to investigate how these associations occur and how endosymbionts may be acquired and transmitted in nature. In our exploratory study of these associations, we found some analytical problems that may be specific for the system, including species that seem more recalcitrant to yield data from their microbiota and a pervasive problem with the use of our choice of V3–V4 primers in an alder-centric environment. However, this experience also helps pointing the way as well as new questions for subsequent approaches to investigate this system, including the use of methodologies that effectively avoid the amplification of organelle sequences (e.g., Hanshew et al. 2013, Beckers et al. 2016), and also investigating causes that could explain the nonrandom efficiency of library preparation comparing Calligrapha with the other species analyzed. We found high and low frequency associations of Wolbachia and Rickettsia with each species investigated, which have been interpreted as stable and transient associations, respectively, with the host. Moreover, we proposed the low frequency, transient associations to be probable because all the species share the same host plant with the particular donor of this strain, Ch. mainensis in this case. However, some elements of this hypothesis need to be verified with population level analyses of these hosts and their endosymbiont associations, also with genetic markers that allow for higher discrimination of bacteria strains, and in localities where these species do not coexist, to exclude the possibility of contemporary, transient transmission. Moreover, as hinted by the results on Ca. philadelphica, other congeneric species with the same and different host plant choices and sympatric with the species studied here would likely provide a useful comparative context for the observed associations with bacteria, including endosymbiont species. Indeed, there are close relatives of A. ambiens and Ch. mainensis that also feed on Alnus (e.g., A. alni Harris, Ch. interrupta F.; Clark et al. 2004), and others that are specialized on Salix (e.g., A. subplicata LeConte, A. bimarginata Say, Ch. knabi Brown; Clark et al. 2004). Relative to the low frequency associations, retrievable thanks to the use of NGS approaches, they have been important for a deeper understanding of the system. In our case, we scored them preferentially as transient stages (possibly in the gut) upon ingestion in the case of the beetle samples, or upon smearing from a Ch. mainensis donor in the case of the leaf samples. These transient stages are recognized as a necessary but insufficient condition to allow for horizontal transmission to occur (Chrostek et al. 2017). However, it would be interesting to verify if the presence of these low frequency strains actually represents a historical infection of these hosts and not contemporaneous to this study. Or, even if the infection were contemporary, discovering whether the infection found its way to the germinal line or to somatic tissues only, which could explain both stable coinfections and differential population densities of the endosymbiont (Pietri et al. 2016). The ecological and evolutionary prospects of these studies are enticing and the analysis of more natural systems like the one we put forward in this work can help a better understanding of the mechanisms behind the association and transmission of endosymbionts to new hosts.