Microbiome analyses of 12 psyllid species of the family Psyllidae identified various bacteria including Fukatsuia and Serratia symbiotica, known as secondary symbionts of aphids

Psyllids (Hemiptera: Psylloidea) comprise a group of plant sap-sucking insects that includes important agricultural pests. They have close associations not only with plant pathogens, but also with various microbes, including obligate mutualists and facultative symbionts. Recent studies are revealing that interactions among such bacterial populations are important for psyllid biology and host plant pathology. In the present study, to obtain further insight into the ecological and evolutionary behaviors of bacteria in Psylloidea, we analyzed the microbiomes of 12 psyllid species belonging to the family Psyllidae (11 from Psyllinae and one from Macrocorsinae), using high-throughput amplicon sequencing of the 16S rRNA gene. The analysis showed that all 12 psyllids have the primary symbiont, Candidatus Carsonella ruddii (Gammaproteobacteria: Oceanospirillales), and at least one secondary symbiont. The majority of the secondary symbionts were gammaproteobacteria, especially those of the family Enterobacteriaceae (order: Enterobacteriales). Among them, symbionts belonging to “endosymbionts3”, which is a genus-level monophyletic group assigned by the SILVA rRNA database, were the most prevalent and were found in 9 of 11 Psyllinae species. Ca. Fukatsuia symbiotica and Serratia symbiotica, which were recognized only as secondary symbionts of aphids, were also identified. In addition to other Enterobacteriaceae bacteria, including Arsenophonus, Sodalis, and “endosymbionts2”, which is another genus-level clade, Pseudomonas (Pseudomonadales: Pseudomonadaceae) and Diplorickettsia (Diplorickettsiales: Diplorickettsiaceae) were identified. Regarding Alphaproteobacteria, the potential plant pathogen Ca. Liberibacter europaeus (Rhizobiales: Rhizobiaceae) was detected for the first time in Anomoneura mori (Psyllinae), a mulberry pest. Wolbachia (Rickettsiales: Anaplasmataceae) and Rickettsia (Rickettsiales: Rickettsiaceae), plausible host reproduction manipulators that are potential tools to control pest insects, were also detected. The present study identified various bacterial symbionts including previously unexpected lineages in psyllids, suggesting considerable interspecific transfer of arthropod symbionts. The findings provide deeper insights into the evolution of interactions among insects, bacteria, and plants, which may be exploited to facilitate the control of pest psyllids in the future.


All 12 Psyllidae species have Carsonella and at least one other symbiont
MiSeq sequencing of the amplicon libraries yielded 46,568-73,470 pairs of forward and reverse reads for the 12 psyllid species (Supplementary Table 1). Denoising and joining of the paired-end reads along with removal of low-quality or chimeric reads resulted in 37,901-63,866 non-chimeric high-quality reads (Supplementary Table 1). Dereplication of these reads resulted in 207 independent sequence variants (SVs), among which only 43 SVs accounted for > 1% of the total reads (Supplementary Table 2). We focused on these 43 SVs, because the targets of the present study were relatively abundant symbionts with close association with the host psyllids, and filtering with the threshold of 1% was shown to be among the most effective and accurate methods to remove potential contaminants derived from environments and experimental reagents [68]. SVs with a relative abundance of less than 1% are collectively categorized as 'others' in Fig. 1, which correspond to 0.16 -3.56% reads in total in each psyllid species (Supplementary Table 2). Notably simple bacterial communities like these have been reported for sternorrhynchan insects with bacteriomes, including aphids, whiteflies, and other psyllid species [24,28,30,37,[64][65][66]69]. All the SVs with a relative abundance of greater than 1% were highly similar to the sequences that were reported to be of insect symbionts (see below). Taxonomic classification by QIIME2 (Supplementary Table 2) followed by independent BLAST searches and phylogenetic analyses showed that all the 12 psyllid species possess distinct lineages of Carsonella (Fig. 1). Because Carsonella has been repeatedly shown to be cospecified with host psyllids [16,19,25,33,34], the phylogenetic relationship of Carsonella is assumed to be useful to infer that of the host psyllids. In the maximum likelihood (ML) tree, the Carsonella sequences from Psyllinae species formed a clade with those of psyllids belonging to the subfamily Ciriacreminae, and the sequence from Epiacizzia kuwayamai formed an independent clade with those of Aphalaroidinae species (Fig. 2). The exclusion of E. kuwayamai from Psyllinae is consistent with the current classification of this species to the subfamily Macrocorsinae. However, these  (Fig. 2), implying their polyphyly as presumed by Burckhardt et al. [2]. However, this branching pattern also lacked robust statistical support (< 50%). Two types each of Carsonella sequences were detected in Cacopsylla peninsularis and Psylla morimotoi ( Fig. 1, Supplementary Table 2). In C. peninsularis, SV19 and SV29 were 99.8% identical (Supplementary Table 2) and formed a clade supported by a bootstrap value of 71% (Fig. 2). SV32 and SV42 from P. morimotoi were also 99.8% identical (Supplementary Table 2), forming a clade supported by a bootstrap value of 97% (Fig. 2). These may reflect sequence variations in each lineage of Carsonella. Although we cannot exclude the possibility that they are artifacts due to polymerase chain reaction (PCR)/sequencing errors, the latter seems less likely because the dada2 plugin corrects sequencing errors during the denoising process [70,71]. Some previous studies that analyzed psyllid microbiomes using 'universal primers' detected only a trace amount of Carsonella reads [27-29, 64, 66, 67]; however, the present study, which used primers appropriately modified to improve sensitivity to highly AT-biased symbiont genes [21,22,31,35], detected a large percentage of Carsonella reads ( Fig. 1, Supplementary Table 2), which reflects actual populations more precisely [30].
Besides Carsonella, all 12 psyllids analyzed in the present study possessed at least one other symbiont ( Fig. 1).

Various Enterobacteriaceae bacteria reside in Psyllidae
Of the 43 SVs obtained in the present study, 39 corresponded to gammaproteobacteria, among which 22 belonged to the family Enterobacteriaceae (order Enterobacteriales) (Supplementary Table 2). Enterobacteriaceae is a group of bacteria that encompasses an especially large fraction of intimate insect symbionts, including those associated with the bacteriome [36,46]. Enterobacteriaceae bacteria identified in the present study include Arsenophonus, Fukatsuia, Serratia, Sodalis, endosym-bionts2, and endosymbionts3. Among them, the most prevalent was endosymbionts3, which is a genus-level monophyletic group of endosymbionts assigned by the SILVA rRNA database project [72].
Whereas Arsenophonus shows a wide range of associations from parasitic to obligately mutualistic to the host insects [74,75], its ecological role in psyllids is currently unknown. with these sequences in the ML tree (Fig. 3). To our knowledge, this is the first formal report of Fukatsuia detected in psyllids. Fukatsuia has only been recognized as a secondary symbiont of aphids, with a wide variety of reported roles, including pathogen, parasite, defensive symbiont, and obligate nutritional symbiont [76]. The recently revealed culturability of Fukatsuia [76] indicates its ability to survive outside aphids, which would facilitate horizontal transfer to other insects, including psyllids. It would be interesting to assess the prevalence and functional role of Fukatsuia in Psylloidea.

Detection of Serratia symbiotica in psyllids
Six SVs found in C. coccinea corresponded to the sequence of Serratia symbiotica, known as a prevalent secondary symbiont of aphids. Namely, SV35, SV37, SV38, SV39, SV40, and SV41, which accounted for 3.0, 2.0, 2.0, 2.0, 1.8, and 1.8% of the denoised C. coccinea reads, respectively, were 98.8 -99.8% identical to a single sequence of S. symbiotica (AB522706) (Fig. 1, Supplementary Table 2). This reference sequence was derived from various aphid lineages, including Acyrthosiphon pisum, Aphis fabae, Aphis gossypii, Cinara pinikoraiensis, Cinara ponderosae, and Trama caudata (all Aphididae). The SVs and S. symbiotica sequence from aphids formed a robustly supported clade (bootstrap: 97%) in the ML tree (Fig. 3). To our knowledge, this is the first formal report of S. symbiotica or its close relative detected in psyllids, although there was a previous mention with no concrete data [77]. These SVs were 98.4% (SV35 vs SV41) -99.8% (SV38 vs SV40) identical to one another. The similarities both in nucleotide sequences and read frequencies imply that the SVs correspond to multiple copies of the 16S rRNA gene in a single S. symbiotica genome. This is consistent with the fact that genomes of several S. symbiotica strains encode more than a single copy of the 16S rRNA gene [42,43], which contrasts the case of primary symbionts with an extremely streamlined genome encoding only a single copy. Similar to Fukatsuia, the ecological role of S. symbiotica is reported to be widely varied depending on aphid lineages [43,44]; however, its role in psyllids is currently unknown. Further studies are required to assess this aspect. As Pons et al. showed that S. symbiotica can enter plants and cause new infection in aphids, host plants are likely media for intraand interspecific horizontal transmission of this bacterium [77].

Sodalis symbionts and its relative
Sodalis endosymbionts were detected in C. burckhardti and C. kiushuensis (Fig. 1, Supplementary [28] and SV34 (see below), whose branching pattern was moderately supported (bootstrap:75%) (Fig. 3). Sodalis symbionts have been detected in a wide variety of insects and are known to have replaced more ancient predecessor symbionts in weevils (Coleoptera: Curculionoidea) [78] and spittlebugs (Hemiptera: Cercopoidea) [79]. In this context, the distribution of the Sodalis symbiont in Cacopsylla spp. may be of interest. Whereas endosymbionts3 appear dominant (presumably bacteriome-associated) secondary symbionts in Cacopsylla spp., C. kiushuensis additionally has a Sodalis symbiont, and C. burckhardti has only Carsonella and Sodalis. This might imply, though speculative, that replacement of endosymbionts3 by Sodalis is at initial stage in C. kiushuensis, and is completed in C. burckhardti. Regarding SV34, which was derived from 6.1% of denoised P. morimotoi reads, QIIME2 failed to assign a genus-level taxonomy (Fig. 1, Supplementary Table 2). The BLAST best hit of SV34 was Sodalis endosymbiont of the psyllid Cardiaspina maniformis (Aphalaridae: Spondyliaspidinae) (KY428659) [28]. These sequences formed a cluster in the ML tree (Fig. 3). However, this branching pattern was only poorly supported (bootstrap: 31%), and their sequence identity was 93.0%, which was below the generally used arbitrary genus threshold of 94.5 -95% [80,81]. Thus, we refrained from assigning this symbiont to a particular genus.

Putative endosymbionts2 symbiont
QIIME2 assigned SV21, which was derived from 31.1% of denoised E. kuwayamai reads, to endosymbionts2 (Fig. 1,  Supplementary Table 2), another monophyletic group of endosymbionts assigned by SILVA [72]. The BLAST best hit of SV21 was a secondary endosymbiont of Cacopsylla myrthi (AF263559) [17], but the sequence identity was only 90.9%. SV21 branched basally to other Enterobacteriaceae bacteria in the ML tree (Fig. 3). It would be interesting to assess the prevalence of endosymbionts2 in the subfamily Macrocorsinae, in the context of the apparent prevalence of endosymbionts3 among Psyllinae species analyzed in the present study.

Psylla morimotoi has Pseudomonas and Diplorickettsia
Non-Enterobacteriales gammaproteobacteria found in the present study were Carsonella (Oceanospirillales: Halomonadaceae) mentioned above, Pseudomonas (Pseudomonadales: Pseudomonadaceae), and Diplorickettsia (Diplorickettsiales: Diplorickettsiaceae); of these, the latter two were detected from P. morimotoi. QIIME2 assigned SV12 and SV30, which were derived from 51.5 and 13.3% of denoised P. morimotoi reads, respectively (Fig. 1, Supplementary Table 2), to Pseudomonas. SV12 and SV30 shared 98.4% identity. SV12 was 100% identical to the sequences of various Pseudomonas strains, including type strains for Pse. graminis (Y11150) and Pse. rhizosphaerae (CP009533). SV30 was 99.5% identical to the sequence of the type strain of Pse. viridiflava (NR_114482). Although Pseudomonas species have been detected in various insects including psyllids, they are largely believed to be transient associates [28]. In contrast to Enterobacteriaceae bacteria, many of which have intimate and stable mutualistic relationships with insect hosts, known examples of Pseudomonas with such associations (vertically-transmitted endosymbionts present in the host hemocoel or cells) are limited in rove beetles (Coleoptera: Staphylinidae) [82] and the adelgid Adelges tsugae (Hemiptera: Sternorrhyncha: Phylloxeroidea: Adelgidae) [83]. Although SV12 and SV30 were not closely related to these symbionts (Fig. 4), the fact that the majority (64.8%) of reads in P. morimotoi corresponded to Pseudomonas (Fig. 1, Supplementary Table 2) implies that the Pseudomonas symbionts potentially play important roles in this psyllid.
QIIME2 assigned SV36, which was derived from 3.7% of denoised P. morimotoi reads, to Diplorickettsia (Fig. 1, Supplementary (Fig. 5). Diplorickettsia massiliensis was observed in Ixodes ricinus and serum samples of human patients with suspected tick-borne disease, suggesting that this bacterium is a human pathogen [84,85]. Subsequently, Diplorickettsia lineages were unexpectedly found in two plant sap-sucking hemipteran insects, M. sexnotatus collected in Japan [86], and D. cf. continua collected in Corsica [30]. The present study adds another example of Diplorickettsia. These findings imply that Diplorickettsia is actually prevalent in various sap-sucking insects. Although their host plants are not shared among M. sexnotatus (Poaceae and Fabaceae), D. cf. continua (Thymelaeaceae), and P. morimotoi (Rosaceae), it would be worth assessing the possibility that the plants are also infected with Diplorickettsia. Diplorickettsia is closely related to the genus Rickettsiella (Diplorickettsiales: Diplorickettsiaceae) (Fig. 5) comprising intracellular bacteria associated with various arthropods, including insects, arachnids, and isopods [19,28]. Whereas many Rickettsiella spp. are simply pathogenic to arthropods, Ca. Rickettsiella viridis [87] modifies the body color of aphids, potentially affecting the attractiveness of insects to natural enemies [88]. As little is known about the functions of Diplorickettsia in host arthropods, it would be interesting to assess the physiological and ecological effects of Diplorickettsia on psyllids.

First detection of Liberibacter in Anomoneura mori
The analysis detected Ca. Liberibacter europaeus (Alphaproteobacteria: Rhizobiales) for the first time in A. mori, a sericultural pest that feeds on the mulberry plants Morus spp. (Moraceae) (Fig. 1, Supplementary Table 2). Namely, SV33, which was derived from 4.6% of denoised A. mori reads, was 99.8% identical to the sequence of Ca.
It appears that Ca. Liberibacter lineages have evolved in close association with Psylloidea, and all known vectors for all Ca. Liberibacter spp. are psyllids [8, 30, 67, 89-91, 94, 95, 97-99]. In this context, the finding that the fecundity and population growth rates of D. citri harboring Ca. L. asiaticus are increased as compared with uninfected insects [101] is particularly interesting. This observed benefit may be an ecological driver for the close association between Ca. Liberibacter spp. and psyllids. Future studies should focus on assessing the general applicability of this hypothesis to other Ca. Liberibacterpsyllid combinations.
Wolbachia are rickettsial bacteria distributed in a wide variety of arthropods and nematodes [104][105][106], whose strains are currently classified into supergroups A-Q [107]. Whereas supergroups A and B are monophyletic and are the most common supergroups infecting arthropods, supergroups C and D infect nematodes. Supergroups E-Q are found in various hosts, including nematodes, springtails, termites, fleas, aphids, and mites [106]. The molecular phylogenetic analysis placed SV7 and SV43 detected in the present study in the robustly supported clade of Wolbachia supergroup B (bootstrap: 97%) (Fig. 7). The majority of Wolbachia strains manipulate the reproduction of arthropod hosts through cytoplasmic incompatibility, feminization, male killing, and parthenogenesis to increase the prevalence of infected females in host populations [104][105][106]. Due to this ability to boost dissemination, Wolbachia are recognized to be promising agents to control insect pests by affecting their traits or microbiomes, including pathogens therein [108,109]. Because of the high infection rates of Wolbachia in pest psyllids worldwide [18,62,66,[110][111][112][113][114][115], and the suggested interactions between Wolbachia and other symbionts [59][60][61][62]116], the application of Wolbachia to control pest psyllids and/or plant pathogens is anticipated [59,62,111,113,115]. The present study suggests rampant horizontal transmissions of Wolbachia among various insect lineages, including pest psyllids, implying the feasibility of artificial infection and/or removal of Wolbachia in psyllids. Such techniques would facilitate the exploitation of Wolbachia as a tool to control pest psyllids and/or the plant pathogens they transmit.
Rickettsia is closely related to Wolbachia, both of which belong to the order Rickettsiales. Similar to Wolbachia, some Rickettsia lineages cause reproductive disorders in host insects, including male killing and parthenogenesis [117,118]. It would be worthwhile to assess whether Rickettsia endosymbionts manipulate the reproduction of psyllids, which will be potentially useful to exploit Rickettsia as a tool to control pest psyllids and/or plant pathogens.

Conclusions
The present study found Ca. Fukatsuia symbiotica and Serratia symbiotica, which were recognized as aphid secondary symbionts, formally for the first time in Psylloidea. The analysis also found the potential plant pathogen, Ca. Liberibacter europaeus (Rhizobiales: Rhizobiaceae), for the first time in a pest psyllid feeding on the mulberry. Furthermore, Wolbachia and Rickettsia, plausible host reproduction manipulators, were detected among analyzed psyllids. The study also identified Arsenophonus, Sodalis, endosymbionts2, endosymbionts3, Pseudomonas, and Diplorickettsia, a plausible human pathogen. These findings suggest considerable interspecific transfer of arthropod symbionts, providing deeper insights into the evolution of interactions among insects, bacteria, and plants. This may be exploited to facilitate the control of pest psyllids with the aid of future studies to determine the localization and genomic structure of the identified bacteria.

Insects and DNA extraction
Adults of 12 psyllid species belonging to the family Psyllidae were collected from several trees of each host plant species in Japan (Table 1). Although most species were selected from the subfamily Psyllinae, Epiacizzia kuwayamai belongs to the subfamily Macrocorsinae [63]. We included this species because it originally belonged to Psyllinae with the name Psylla kuwayamai and its generic transfer from Psylla to Epiacizzia appeared to be inconsistent with the morphological features [119]. Also, whereas Burckhardt et al. recently proposed to transfer Psylla morimotoi (Psyllinae) to the genus Spanioneura (Psyllinae) based on its ecological features [2], the present study refers to this species as is because its morphological features are more consistent with the definition of Psylla (sensu stricto) [120,121].
After surface sterilization with ethanol, DNA was extracted from the whole bodies of pooled individuals of five adult males and five adult females of each psyllid species using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). The quality of the extracted DNA was assessed using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts, U.S.A.) and the quantity was assessed using a Qubit 2.0 Fluorometer with the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific).

Construction and sequencing of amplicon libraries
Bacterial populations in psyllids were analyzed using the MiSeq system (Illumina, San Diego, California, U.S.A.), as described previously [30]. Briefly, amplicon PCR was performed using the genomic DNA extracted from psyllids, the KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, Massachusetts, U.S.A.), and the primer set 16S_341Fmod (5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG YYTAMGGRNGGC WGC AG-3′) and 16S_805R (5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGA CTACHVGGG TAT CTA ATC C-3′) targeting the V3 and V4 regions of the 16S rRNA gene. Whereas both primers were based on the instructions by Illumina [122], 16S_341F was modified (underlined), where original CC, C, and G were replaced with mixed bases, YY (C or T), M (A or C), and R (A or G). Our preliminary analyses demonstrated that this modification improves sensitivity to detect symbionts with ATrich genomes including Carsonella, without reducing sensitivity for those with GC-rich genomes. Dual indices and Illumina sequencing adapters were attached to the amplicons with index PCR using the Nextera XT Index Kit v2 (Illumina). The libraries were combined with the PhiX Control v3 (Illumina), and 250 bp of both ends were sequenced on the MiSeq platform (Illumina) with the MiSeq Reagent Kit v2 (500 cycles; Illumina).

Computational analysis of bacterial populations
After the amplicon sequence reads were demultiplexed, the output sequences were imported into the QIIME2 platform (version 2020.2) [123] and processed as described previously [30]. In brief, after primer sequences were removed using the cutadapt plugin [124], paired-end sequences were trimmed, denoised, joined, and dereplicated using the dada2 plugin [70]. During this step, chimeric sequences were detected and removed. The q2-feature-classifier plugin [125] was trained using the V3 and V4 regions of the 16S rRNA gene sequences retrieved from the SILVA database ver. 132 (SILVA_132_ QIIME_release/taxonomy/16S_only/99/taxonomy_7_levels. txt) that were clustered at 99% sequence similarity [72]. Subsequently, the denoised and dereplicated amplicon reads were classified and taxonomic information was assigned using the trained q2-feature-classifier. Obtained sequence variants (SVs) were manually checked by performing BLASTN searches against the National Center for Biotechnology Information non-redundant database [126].

Phylogenetic analysis of detected bacteria
Phylogenetic analysis of SVs was performed as described previously [30]. Briefly, after the SVs were aligned to related sequences using SINA (v1.2.11) [127], phylogenetic trees were inferred by the maximum likelihood (ML) method using RAxML (version 8.2.12) [128]. The GTR + Γ model was used with no partitioning of the data matrix, with 1000 bootstrap iterations (options -f a -m GTRGAMMA -# 1000).