Aphid Species Specializing on Milkweed Have Distinct Bacterial Symbiont Communities


 Host plant range is arguably one of the most important factors shaping microbial communities associated with insect herbivores. However, it is unclear whether host plant specialization limits microbial community diversity or to what extent herbivores sharing a common host plant evolve distinct microbiomes. To investigate whether variation in host plant specialization influences the composition of herbivore symbiont populations we compared bacterial diversity across three milkweed aphid species (Aphis nerii, Aphis asclepiadis, Myzocallis asclepiadis) feeding on a common host plant (Asclepias syriaca) using 16S rRNA metabarcoding. Overall, bacterial species richness did not vary with degree of host plant specialization. However, aphid species harbored distinct bacterial communities that varied in composition and relative abundance of key symbionts. Differences in aphid microbiomes were primarily due to strain variation in the obligate symbiont Buchnera and facultative symbiont Arsenophonus, as most of the low-abundant taxa were found in all three species. Interestingly, A. asclepiadis harbored a greater diversity of unique strains of Buchnera and significantly higher Arsenophonus relative abundances compared to the other two aphid species. Although many low abundance microbes were shared across all milkweed aphids, key differences exist in symbiotic partnerships that could influence additional ecological variation, including variation in ant tending observed across milkweed aphid species via microbial induced changes to honeydew or defensive chemical profiles. This study suggests generalist and specialist herbivore microbiomes are similar when feeding on a common host plant and highlights an intriguing potential role for strain level variation of key aphid symbionts in host-plant interactions.


Introduction
Establishing partnerships with microbes is hypothesized to enhance the evolutionary and adaptive potential of the host organism [1][2][3][4]. In fact, both plants and insects harbor diverse microbial communities that provide unique advantages, including protection from stress and novel mechanisms for nutrient acquisition [4][5][6][7]. However, only recently have microbes come to be viewed as the architects of plant-insect interactions. From an insect perspective, microbes can be critical in facilitating or restricting the use of host plants by aiding in digestion of plant material and detoxi cation of anti-herbivore chemical defenses [3,5,8]. Soil and root-associated microbes also play essential roles in plant growth and defense against both above and belowground insect attackers [6,7,9,10]. Although microbes unquestionably impact plant-insect interactions, the evolutionary forces shaping plant and insectassociated microbial communities remain poorly understood.
One factor thought to play a pivotal role in shaping insect herbivore microbiomes is host plant range or diet breadth, which can vary from a single species to hundreds of plant species [11][12][13][14][15][16]. Several plant associated factors, including nutritional quality and defensive chemistry, are known to in uence insect microbial community dynamics by selecting for taxa that facilitate enhanced colonization and survival [17]. For example, associations with heritable bacterial symbionts have long been recognized as enabling phloem feeding insects to exploit an otherwise poor nutritional source [17,18]. Microbes are also hypothesized to contribute to variation in the capacity of insect herbivores to consume chemically defended plants [8]. Finally, feeding on different host plants can in uence the composition of microbial communities associated with insect herbivores [11,14,19,20]. However, aside from a handful of wellcharacterized nutritional symbionts it is unclear to what extent insect microbial communities either directly or indirectly contribute to host plant adaptation and diet breadth. Furthermore, the mechanisms contributing to variation in herbivore microbiomes and ecological implications for plant-insect coevolution are not well understood [14].
Insect microbiomes consist of varying combinations of heritable symbionts and exible pools of environmentally acquired microbes, but how host plant range shapes the taxonomic and functional diversity of these communities is relatively unexplored [5,[21][22][23]. Microbial communities of specialist herbivores may be less diverse relative generalist species that possibly require a larger microbial repertoire to successfully colonize highly variable hosts or simply acquire greater environmental microbial diversity from feeding on a broader host plant range [5,21]. Alternatively, overlapping host plant range could result in selection for common microbes needed for successful colonization across different herbivore species due to similarities in nutritional ecology or defensive chemistry. Currently, it is unknown whether broader patterns associated with host plant use exist, such as a gradient in microbiome diversity from specialist to generalist herbivores.
The current study aims to address the following question: How does the degree of host plant specialization shape insect herbivore symbiont community diversity and composition? Heritable bacterial symbionts are hypothesized to expand the diet breadth of sap-feeding insects, particularly aphids [23][24][25]. However, aphid studies have primarily focused on the role microbial symbionts play in host plant use by highly polyphagous species, while much less is known about monophagous or oligophagous species.
Aphids colonizing milkweed species in the subfamily Asclepiadaceae show varying degrees of host specialization [26] and thus represent an excellent system in which to investigate the role diet breadth plays in shaping aphid microbial community composition. Aphis nerii is considered a broad specialist feeding on over 50 plant species, including oleander, milkweed, and periwinkle. Aphis asclepiadis is a narrow specialist, feeding on less than 10 Asclepias spp., and Myzocalis asclepiadis is monophagous, feeding only on the common milkweed Asclepiassyriaca [26,27]. It is thought that milkweed aphid species have few microbial symbionts besides the universal obligate symbiont Buchnera, possibly due to the toxicity of milkweed defensive chemicals (e.g. cardenolides) which when ingested could make the internal aphid environment inhospitable for sustained microbial growth (e.g. via direct anti-microbial effects). However, few populations have been screened for symbionts and thus the bacterial community of milkweed aphids remains largely unexplored [28]. We therefore characterized the diversity and composition of bacterial communities using targeted amplicon sequencing of eld collected aphids found naturally colonizing common milkweed (A. syriaca).

Aphid Field Collection & Sample Preparation
Aphids were collected in July-Aug 2017 from a single host plant type, the common milkweed (A. syriaca), across 14 locations generally surrounding the Purdue campus (Online Resource 1 Table S1). Locations were chosen based on the presence of large patches of common milkweed in open grassy elds away from roadsides. Distance between locations ranged from a minimum of 1km to a maximum of 141km (Online Resource 1 Table S1). To avoid sampling heavily from single aphid clonal families (e.g. offspring from a single alate aphid colonizing a plant), aphids were collected from 3-5 separate plants (<20m apart) within a single location (e.g. from distinct milkweed patches when possible) and pooled together into a single sample. In some cases, 2 aphid species were found co-colonizing a single plant and thus single species versus mixed species samples were designated as such during data collection (Online Resource 1 Table S1). Individual aphids were removed from milkweed leaves by hand using a paint brush, identi ed and separated by species, and then placed in Eppendorf storage tubes containing 95% ethanol. All individuals collected from each separate plant at each location were stored at -20C until further processing. The total number of aphids collected per species at a single location ranged from 10 individuals to >50 individuals.
For bacterial community pro ling, groups of 5 individual aphids per species were selected from the total pooled sample per location for DNA extraction and sequencing. Total DNA was extracted from groups of 5 whole aphids using the Qiagen DNeasy kit following standard protocols. Aphids were not surface sterilized beyond storage in 95% ethanol and therefore, while bacterial DNA recovered is expected to be primarily from internal communities, not all external or cuticular bacterial DNA was eliminated and is expected to contribute to overall diversity estimates [29]. Total DNA concentration was measured using a Nanodrop spectrophotometer for all samples. In total 65 pools of aphids were processed and sequenced for microbiome analysis. Targeted amplicon sequencing was used to characterize the bacterial communities associated with the three milkweed aphid species (A. nerii, A. asclepiadis, M. asclepiadis) across the 19 locations sampled. Library preparation and sequencing was performed at the University of Minnesota Genomics Core Facility on an Illumina MiSeq instrument (V3 cluster chemistry) following optimized methods [30] and using primers that target a 250bp segment of the V4 region of the 16s rDNA gene. Primers used were standard V4 region primers 515F-GTGCCAGCMGCCGCGGTAA and 806R-GGACTACHVGGGTWTCTAAT. All sequence data is available on NCBI SRA database under project number PRJNA635683. and low quality reads. All subsequent processing was performed in R (v 3.6.3) and Bioconductor (v 3.10). Trimmed reads were processed through the dada2 (v 1.14.1) pipeline by ltering and trimming based on read quality, inferring error rates, merging paired end reads, removing chimeras, and assigning taxonomy with the Silva reference database v. 132 [33]. Removal of very low abundance reads was done using a cut-off of fewer than 10 reads in 5% of the samples. Next, eukaryotic and mitochondrial sequences were removed. Lastly, samples with fewer than 2500 reads were removed. A summary of raw sequencing results and processing steps can be found in Online Resource 1.

Characterization & Analysis of Aphid Bacterial Communities
Following sequence processing, all downstream analyses were run in R (v 3.6.3). All code for statistical analyses and generation of gures can be found in the Purdue University Github (https://github.itap.purdue.edu/LaramyEndersGroup/Milkweed-Aphid-Microbiome). Speci cally, we compared standard alpha and beta diversity metrics using the phyloseq [v 1.30.0; 32] and vegan [v 2.5-6; 33] packages in R to determine the extent to which microbiomes varied in taxonomic composition and structure across aphid species and sampling locations. To compare species richness and evenness we used the Kruskal-Wallis test to test the effect of aphid species on diversity metrics, followed by the Wilcoxon Rank Sum test to identify signi cant differences among aphid species. Differences in the structure of bacterial communities across aphid species was assessed through PERMANOVA analysis of beta diversity (Unifrac, weighted Unifrac, Bray-Curtis) and visualized using Non-Metric Multidimensional Scaling (NMDS). Homogeneity of dispersion across aphid species' samples was also tested using PERMDISP. To identify differentially abundant bacterial sequences across aphid species we applied generalized linear mixed models to normalized read counts using DeSeq2 [v 1.26.0; 34] and made pairwise comparisons between each aphid species. Normalized read counts account for variation in library size and were used to estimate relative abundance of each ASV. Finally, we compared strain level genetic differences in Buchnera and Arsenophonus ASVs using a phylogenetic approach by aligning sequences with the DECIPHER package [v 2.14.0; [37] and building Maximum-likelihood trees with the optim.pml(model="GTR", rearrangement = "stochastic") function in the phangorn package in R [v 2.5.5; 35].

Results
Overall, most bacterial taxa identi ed through sequencing were found to occur in all three milkweed aphid species (Fig 1a). In total 49 ampli ed sequence variants (ASVs) were identi ed (Online Resource 1 Table  S3), 42 of which occurred in each species when all sampling locations were considered. Bacterial communities in general were dominated by the primary aphid symbiont Buchnera, which was in high relative abundance compared to the remaining taxa (Fig 1b). Among the non-Buchnera ASVs identi ed several well-known aphid facultative symbionts were found, including Arsenophonus, Serratia, and Hamiltonella, which occurred in all milkweed aphid species to varying degrees (Fig 1b). In a few samples Regiella, Ricketsiella, and Wolbachia were identi ed in very low abundance (<2 reads in 10% of samples) and did not make the cutoff to be included in analysis (Online Resource 2 Fig S1). Much of the remaining low-abundance taxa identi ed are likely environmentally acquired from the plant host or surrounding soil (e.g. Pantoea, Pseudomonas, Streptococcus). The prevalence or proportion of samples that tested positive for each ASV also tended to vary across aphid species (Online Resource 1 Table S3). For example, Buchnera was found in all samples, but Hamiltonella ranged from 13% (2/15) prevalence in A. asclepiadis samples to 100% (6/6) in M. asclepiadis and 30% (13/43) in A. nerii samples. Serratia was also highly prevalent across all three species, ranging from 83% in A. nerii to 100% in M. asclepiadis samples. Interestingly, all three identi ed strains of Lactobacillus were found in 100% of M. asclepiadis samples but were only present in 50% or less of samples of the other two species. Co-occurring strains of Buchnera and Arsenophous were also found in all three species, although it should be noted that we tested groups of aphids and it remains unknown whether individual aphids harbor different strains (i.e. multiple infections).
Although most ASVs were shared at a global level, several key differences indicate that each milkweed aphid species harbors distinct bacterial communities that vary in composition and abundance (Figs 1 & 2, Online Resource 1 Table S4). NMDS of Bray-Curtis dissimilarity between samples indicates that bacterial community structure is unique to each aphid species (p <0.001, Online Resource 1 Table S5), but did not vary by location (Fig 2 & 3). Hierarchical clustering of samples by differences in community structure and relative abundance of each taxa further shows that each aphid species forms a unique group (Fig 3). In addition to differences in overall community structure there was signi cant variation in the bacterial species richness and evenness of communities found across the three aphid species (Fig 4, Online Resource 1 Table S5). Interestingly, A. nerii had the lowest bacterial species richness and A. asclepiadis had the highest among the three species. Overall, results did not show a gradient in bacterial diversity (i.e. species richness) associated with diet breadth or degree of host-plant specialization across these three milkweed aphid species.
To further identify factors contributing to differences between aphid microbiomes we compared the relative abundance of individual taxa and variation in bacterial strain diversity (ASVs) (Online Resource 1 Table S5; Figures 3 & 5). Differences in community structure were primarily driven by 1) variation in the composition of facultative symbionts and other low-abundance (non-Buchnera) bacteria (Figs 1b & 3) and 2) the presence of Buchnera and Arsenophonus strains unique to each aphid species (Fig 5). Overall, A. nerii bacterial communities were dominated by Buchnera with few other symbionts in low abundance, which differs from the other two aphid species. A. asclepiadis has a greater diversity of unique Buchnera strains and overall higher Arsenophonus abundances. Each aphid species typically harbored a single dominant Buchnera strain in highest relative abundance, with some co-occurring strains at lower abundances (Figs 4 & 5; Online Resource 1 Table S3). Where communities differed the most was in Arsenophonus strain diversity and abundance (Figs 3 & 5). For example, there were 6 distinct Arsenophonus strains found within A. asclepiadis communities that ranged between 50 -227 fold higher relative abundance compared to A. nerii and M. asclepiadis. A. asclepiadis also harbored 9 distinct Buchnera strains across the locations sampled. In addition, Buchnera strains varied from 91.6% sequence similarity to several strains differing by only 1-2 base pairs (>99.2% similarity) (Online Resource 1; Table  S7a). Arsenophonus strains were even more similar to each other, varying from 94.7 -99.6% similarity (Online Resource 1; Table S7b).

Discussion
The prevailing view of host-associated microorganisms and their role in insect-plant coevolution is rapidly changing. Acquisition of bene cial microbes serves not only to expand the ecological niche of the host, but can add novel weaponry in the adaptive battle between plants and the herbivores that prey upon them [3,5,7,8]. Furthermore, insects not only acquire microbes from host plants and the surrounding environment, but these microbial communities are intimately linked [9,39,40]. Disentangling the mechanisms driving variation in herbivore microbial communities and the ecological consequences for host plant specialization is therefore of great interest. In the current study, shifts were observed in the core microbiome of aphid species exhibiting varying degrees of host specialization within the milkweed family. Aphids feeding on milkweed share most bacterial taxa (Fig 1), but each species harbors unique symbiont populations that differ in strain diversity and relative abundance (Figs 1b, 3 & 5). However, broader patterns in bacterial symbiont diversity that scale with diet breadth were not observed (Figs 3 &  4).
Host plant speci city can in uence herbivore associated microbial communities on multiple levels, including causing changes in functional, taxonomic or strain level diversity or by altering the abundance of individual taxa. The microbiomes of generalists and specialist herbivores could also vary simply due to differences in the contribution of heritable versus environmentally acquired microbes, the latter being more variable and transient. However, diet may in fact generally be a poor predictor of insect bacterial community composition [41], with exceptions only in certain groups of herbivores [42]. For example, Lepidopteran larva have variable gut microbiomes that are likely shaped by the host plant they feed on, but general patterns associated with diet breadth are not observed due to the high turnover rate of these microbial communities [14]. In contrast, recent work in Costa Rican rolled-leaf beetles (Cephaloleia spp.) shows diet breadth is linked to microbiome diversity and community structure [21]. Previous research also indicates the composition of aphid microbial communities are structured by host plant [11,15,16,20,22,24,43,44] and heritable bacterial symbionts are involved in expanding diet breadth [23,25].
Our results show milkweed aphids differing in host plant specialization and diet breadth have similar microbial communities on a taxonomic level but differ in strain diversity and relative abundance of key symbionts (e.g. Arsenophonus). One possible explanation is exposure to similar nutritional and chemical pro les could homogenize microbiomes of herbivore species feeding on the same host plants. The species in this study are specialists of the Asclepias (milkweed) family and therefore may have similar core symbionts due to exposure to closely related host plants. Additionally, phylogenetic relatedness can generally result in closely related species harboring more similar microbiota than distantly related species [45]. The current study is limited in that aphids were collected from a single host plant (i.e. common milkweed) and generalist species with a host range outside the milkweed family (e.g. Myzus persicae) were not characterized due to low occurrence in the eld. Additional imbalance in sampling across aphid species (e.g. fewer M. asclepiadis samples, Table S1) could mean that some microbial variation was missed. Further research investigating the generalist-specialist gradient using herbivores that feed across multiple plant families is needed to clarify the extent to which diet breadth shapes microbial communities (e.g. [21]).
Although milkweed aphid microbiomes were overall similar, key differences in symbiotic partnerships could contribute to additional ecological variation (e.g. ant tending, parasitism rates, predation). Interestingly, variation in symbiont relative abundance and strain diversity contributed most to differences observed across milkweed aphid microbiomes. In particular, Arsenophonus was found in higher abundance in A. asclepiadis compared to the other two species. Arsenophonus is a notorious shapeshifting insect symbiont, known best for reproductive manipulation of its host [46,47]. Most aphid facultative symbionts are found in much lower abundances compared to the obligate nutritional symbiont Buchnera, suggesting the unusually high abundance observed in A. asclepiadis could be linked to symbiont complementarity, as has occurred in other aphid species [48]. It is also possible Arsenophonus provides a general tness boost, similar to what has been observed in the soybean aphid [49]. Finally, differences in symbiont populations could shape milkweed aphid-ant mutualisms, possibly via microbial induced changes to honeydew or emission of chemical compounds that mediate partner attraction. Previous work shows insect social partnerships not only uniquely in uence each host's symbiotic microbiome [50], but that volatile organic compounds produced by aphid-associated microbes play a role in attracting ant mutualists [51]. Among milkweed aphids, A. asclepiadis is consistently tended by ants and bene ts from enhanced protection from predators, while A. nerii is occasionally ant tended and M. asclepiadis appears to be a loner lacking ant friends [52]. Based on the current study an intriguing question arises; Does Arsenophonus mediate milkweed aphid-ant interactions and thus contribute to observed differences in ant attendance? Although Arsenophonus does not appear to in uence the intensity of ant attendance in cowpea aphids [53], it is possible this symbiont has evolved a different ecological role in the case of A. asclepdiadis. However, given the 16s rRNA metabarcoding approach used in this study only provides relative abundances, Arsenophonus titer levels will need to be con rmed using additional methods such as quantitative PCR in order to take the rst step towards addressing a potential role in aphid-ant interactions.
In summary, we did not nd evidence for a gradient in bacterial diversity across generalist and specialist milkweed aphids. Instead, our results suggest overlapping host plant range and shared hosts can result in selection for common microbes and thus highly similar microbiomes across species. However, milkweed aphids do harbor unique bacterial populations that most notably vary in strain diversity and relative abundance of Arsenophonus, although a handful of other well known aphid symbionts were also found in low abundance. These ndings suggest that while host plant specialization may not be a major driver of divergence in overall taxonomic composition of aphid symbiont communities, factors such as strain level variation and differences in abundance offer alternative routes to generating adaptive potential. Further research is needed to determine the functional or ecological role played by milkweed aphid facultative symbionts and different co-occurring strains. Figure 1 A) Venn Diagram showing the extent to which amplicon sequence variants (ASVs) are shared and which are unique to each milkweed aphid species from samples collected across all locations. B) log10(x+1) normalized counts of bacterial taxa found across the three milkweed aphid species. If multiple strains were present they were grouped together under one genus. Very low abundant ASV (below 1% relative abundance) were grouped together into the "Other" category, except for known aphid symbionts (i.e. Hamiltonella) and signi cantly differentially abundant taxa (see Online Resource 1 & 2; Table S3 &   Genetic variation and differences in relative abundance of Buchnera and Arsenophonus strains (ASVs) identi ed in the three milkweed aphid species. All Buchnera and Arsenophonus ASVs except ASV 49 were differentially abundant (adjusted p < 0.05) in at least one species comparison by DESeq2 (Online Resource 1