Aphid species specializing on milkweed harbor taxonomically similar bacterial communities that differ in richness and relative abundance of core symbionts

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 similar microbiomes. To investigate whether variation in host plant range influences the assembly of core 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, although there was significant overlap in taxa detected across all three aphid species (i.e. similar composition), some structural differences were identified within communities. Each aphid species harbored bacterial communities that varied in terms of richness and relative abundance of key symbionts. However, bacterial community diversity did not vary with degree of aphid host plant specialization. Interestingly, the narrow specialist A. asclepiadis harbored significantly higher relative abundances of the facultative symbiont Arsenophonus compared to the other two aphid species. Although many low abundance microbes were shared across all milkweed aphids, key differences in symbiotic partnerships were observed that could influence host physiology or additional ecological variation in traits that are microbially-mediated. Overall, this study suggests overlap in host plant range can select for taxonomically similar microbiomes across herbivore species, but variation in core aphid symbionts within these communities may still occur.

a handful of well-characterized 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 flexible pools of environmentally acquired microbes [21][22][23] , but how host plant range shapes the taxonomic and functional diversity of these communities is relatively unexplored 5,[24][25][26] . Variation in host plant range can influence the assembly of microbial communities and dynamics of symbiotic partnerships via several hypothesized mechanisms. For example, microbial communities of specialist herbivores are hypothesized to be less diverse relative to generalist species that possibly require a larger microbial repertoire to successfully colonize highly variable hosts (e.g. adaptive advantage to harboring multiple heritable symbionts) 5,24 . Generalists may also be exposed to and therefore acquire greater environmental microbial diversity from feeding on a broader host plant range. Alternatively, overlapping host plant range across herbivore species could result in selection for common microbes needed for successful colonization due to similarities in nutritional ecology or defensive chemistry of shared host plants (i.e. purifying selection reduces differences in symbiont communities). Currently, it is unclear whether broader patterns associated with host plant use exist, such as a gradient in microbiome diversity from specialist to generalist herbivores or shifts in abundance of key heritable symbionts.
Heritable bacterial symbionts are hypothesized to expand the host plant range of sap-feeding insects, particularly aphids [26][27][28] . Signatures of host plant specificity have been detected in aphid microbiome composition 20,[29][30][31][32] . Even within a single aphid species ecologically divergent biotypes specialized on different host plants can harbor distinct microbiomes that differ in taxonomic composition and frequency of symbionts 25,33,34 . In contrast, recent work suggests aphid phylogeny (i.e. species relatedness) and geographic isolation of populations are dominant factors predicting differences in bacterial communities, while host plant range may be less important overall [21][22][23] . 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. As a result, the relative importance of ecological factors such as host plant range and geographic range in governing aphid symbiont community assemblages remains unclear.
The current study aims to address the following question: How does variation in host plant range shape the assembly of aphid symbiont communities? Aphids colonizing milkweed species in the subfamily Asclepiadaceae exhibit variation in host plant specialization 35 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 capable of 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 Asclepias syriaca 35,36 . Currently, few populations have been screened for symbionts and thus the bacterial community of milkweed aphids remains largely unexplored 37 . One possibility is that milkweed aphid species have few microbial symbionts besides the universal obligate symbiont Buchnera 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, it is generally unclear how ecological factors (e.g. host plant range, geographic location) shape both the heritable and environmentally acquired microbiota of milkweed aphids. We therefore characterized the diversity and composition of bacterial communities using targeted amplicon sequencing of field collected aphids found naturally colonizing common milkweed (A. syriaca). Specifically, we tested two alternative hypotheses regarding the role host plant range plays in shaping aphid symbiont community assembly. First, symbiont diversity (i.e. bacterial species richness) is hypothesized to increase with expansion of host plant range or increased diet breadth (A.nerii > A. asclepiadis > M. asclepiadis). Alternatively, overlap in host plant range is hypothesized to reduce differences in symbiont communities across aphid species (e.g. via purifying selection from exposure to similar plant defensive chemicals). Overall, we find that milkweed aphids collected from a common host plant tend to have highly similar bacterial microbiomes when considering the broader pool of taxa detected across populations, but differences in relative abundances and strain diversity were found in core heritable symbionts.

Methods
Aphid field collection & sample preparation. In this study we focus on three aphid species (A.nerii, A. asclepiadis, M. asclepiadis) that colonize plants of the family Apocynaceae with varying degrees of host specialization 35,36 . The broad specialist species A.nerii is considered an obligate parthenogen (only reproduces asexually in the wild) and in its introduced range in the United States (U.S.) it colonizes species of Apocynoideae (oleander) and Asclepiadaceae (milkweed) 38,39 . The wide host plant distributions for A. nerii overlap and extend into northern and midwestern regions of the U.S. 40,41 . Therefore populations of A.nerii can potentially feed on multiple host plant species throughout the summer in the midwestern and eastern U.S. However, A. nerii cannot tolerate freezing temperatures and is unable to support year round populations except in the southernmost regions of the U.S. 38,39 . The narrow specialist A. asclepiadis is also broadly distributed across the eastern and central U.S., but it is unknown if this species is an obligate or cyclical parthenogen or whether it can overwinter in the egg stage in northern regions 39 . Even less is known about the life history of the monophagous species M. asclepiadis, which may be cyclically parthenogenic and overwinter in the egg stage. Interestingly, M. asclepiadis is less gregarious compared to the other two species and highly mobile due to the asexual production of only winged adults during the summer 36,42 .
To address our central question and test hypotheses regarding differences in aphid microbiomes associated with degree of host plant specialization, we sampled aphids from a single common host plant across different locations. By sampling from a single host plant type we focus on detecting differences in heritable symbionts www.nature.com/scientificreports/ associated with variation in aphid diet breadth, and to a lesser extent environmentally acquired microbes associated with different locations or individual host plants. Our sampling design avoids confounding effects associated with variation in microbiota that could be acquired from feeding on different host plant species, and instead focuses on detection of variation in the heritable fraction of the aphid microbiome (e.g. broader host plant ranges may select for greater variation in heritable facultative symbiont populations). Further, M. asclepiadis is monophogus and it is not possible in this system to sample all three aphid species from multiple shared host plant species, an experimental design that may be feasible in other systems to further tease apart host plantaphid-microbiome interactions. 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 (Supplementary Table S1) in West Lafayette, IN. Ascelpiadis syriaca is among the most abundant milkweed species in the midwestern United States 43 . Locations were chosen based on the presence of large patches of common milkweed in open grassy fields away from roadsides. Distance between locations ranged from a minimum of 1 km to a maximum of 141 km (Supplementary  Table S1). To avoid sampling heavily from single aphid clonal families (e.g. offspring from a single alate aphid colonizing a plant) that could limit the ability to capture natural variation in microbial communities, aphids were collected from 3 to 5 separate plants (< 20 m apart) within a single location (e.g. from distinct milkweed patches when possible) and pooled together into a single sample for downstream microbiome analysis. 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 (Supplementary Table S1). Individual adult aphids were removed from milkweed leaves by hand using a paint brush, identified and separated by species using standard identification guides [ 35 , http:// www. aphid sonwo rldsp lants. info/], and then placed in Eppendorf storage tubes containing 95% ethanol. Only adult aphids were included because they have distinct morphological characteristics that allow for straightforward species identification (e.g. body coloration; M. acelepiadas has distinctive red-orange blotches and black markings), thus preventing mixed species samples and ensuring pools of aphids used for downstream analysis were separate species. A.nerii and A. asclepiadis reproduce parthenogenetically during the summer and can produce winged individuals under crowded conditions 39 . Therefore, only wingless apterous adults were collected from A. nerii and A. asclepiadis, which are relatively sessile and typically remain on a single plant from birth to death unless disturbed. In contrast, M. asclepiadis only produce alates or winged individuals in the summer, which can move between individual plants of their only host A. syriaca 36 . Overall, individuals used in this study are unlikely to have feed on multiple host plant species prior to field collection. All individuals collected from each separate plant at each location were stored in 95% ethanol at −20C until further processing. The total number of individuals collected per species at a single location ranged from 10 individuals to > 50 individuals.
For bacterial community profiling, 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 surface sterilized with 95% ethanol (during storage) and washed with ultra pure water prior to extraction. While bacterial DNA recovered is expected to be primarily from internal communities, it is possible not all external or cuticular bacterial DNA was eliminated and is expected to contribute to overall diversity estimates 44 . Total DNA concentration was measured using a Nanodrop spectrophotometer for all samples. In total 82 pooled aphid samples were processed and sequenced for microbiome analysis (A. nerii n = 46, A. asclepiadis n = 30, M. asclepiadis n = 6; Supplemental Table S1). 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, paired end 250 bp sequencing) following optimized methods that target the V4 region of the 16 s rDNA gene 45 . Primers used were standard V4 region primers 515F-GTG CCA GCMGCC GCG GTAA and 806R-GGA CTA CHVGGG TWT CTAAT 46 48 ] were used to remove adapters and primer sequences 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; 49 ] pipeline by filtering 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. 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, individual samples with fewer than 2500 reads were removed and a final sample-level filtering step (i.e. Buchnera ASVs with < 1% reads within a sample were removed from individual samples) was applied to the remaining 64 samples to account for potential read contamination or "cross-talk" among samples that can occur in metabarcoding studies of microbial communities that are dominated by few high abundance symbionts (e.g. Buchnera-50 ). We applied the more conservative exclusion cut-off to Buchnera ASVs only in order to avoid loss of rare community members and minimize false positives resulting from mis-binning Buchnera sequences. A summary of raw sequencing results and processing steps can be found in Supplementary Table S2.
Following sequence processing, all downstream analyses and data visualization was run in R (v 3.6.3). All code for statistical analyses and generation of figures, including information on R packages used, can be found in the Purdue University Github (https:// github. itap. purdue. edu/ Laram yEnde rsGro up/ Milkw eed-Aphid-Micro biome). Specifically, we compared standard alpha and beta diversity metrics using the phyloseq [v 1. 30 www.nature.com/scientificreports/ and structure across aphid species and sampling locations. Initial data assessment (i.e. Q-Q plot of residuals, Shapiro-Wilk Normality test and Levene test for homogeneity of variance) indicated our alpha diversity data did not fit the assumptions of a linear model. To compare species richness and evenness we therefore used the non-parametric Kruskal-Wallis test to test the effect of aphid species on diversity metrics, followed by the Wilcoxon Rank Sum test with Holm correction for multiple testing to identify significant 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; 53 ] and made pair-wise 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; 54 ] and building Maximum-likelihood trees with the optim.pml(model = "GTR", rearrangement = "stochastic") function in the phangorn package in R [v 2.5.5; 55 ].

Results
Overall, most bacterial taxa identified through sequencing were found to occur in all three milkweed aphid species (Fig. 1a). In total 45 amplified sequence variants (ASVs) were identified (Supplementary Table S3), 38 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 identified 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 identified in very low abundance (< 2 reads in 10% of samples) and did not make the cutoff to be included in analysis (Supplementary Fig. S1). Much of the remaining low-abundance taxa identified 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 (Supplementary 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 identified 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 Arsenophonus were also found to varying degrees 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 when the entire pool of microbes was considered across locations (Fig. 1a), several key differences indicate that each milkweed aphid species harbors distinct bacterial communities that vary in structure and abundance of taxa (Figs. 1b, 2, Supplementary Table S4). NMDS of   Table S5). 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 significant variation in the bacterial species richness and evenness of communities found across the three aphid species (Fig. 4, Supplementary Table S5). Interestingly, the broad specialist A. nerii had the lowest bacterial species richness and the more narrow specialist 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) (Supplementary Table S6; Figs. 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. Unique to A. asclepiadis was the overall higher Arsenophonus abundances, but lower Buchnera abundances compared to the other two species (Fig. 1b). 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; Supplementary  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 and 227 fold higher relative abundance compared to A. nerii and M. asclepiadis. A. asclepiadis also harbored 5 distinct Buchnera strains across the locations sampled. In addition, Buchnera strains varied from 92.1% sequence similarity to several strains differing by only 1-2 base pairs (> 99.2% similarity) (Supplementary Table S7). Arsenophonus strains were even more similar to each other, varying from 94.7 to 99.6% similarity (Supplementary Table S7).

Discussion
The prevailing view of host-associated microorganisms and their role in insect-plant coevolution is rapidly changing. Acquisition of beneficial 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 colonize 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,56,57 . 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 when the broader pool of microbes detected across populations is considered (Fig. 1), but each species harbors unique populations that differ in strain diversity and relative abundance of core heritable symbionts, most notably  (Figs. 1b, 3, 5). However, broader patterns in bacterial symbiont diversity that scale with diet breadth were not observed (Figs. 3, 4). Instead, our results provide support for the hypothesized effect of overlapping host plant range and common selective pressures (e.g. plant chemical defenses) leading to similarities in microbiomes across milkweed aphid species. Host plant specificity can influence 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 22,58 , with exceptions only in certain groups of herbivores 59 . 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 24 . Previous research also indicates the composition of aphid microbial communities are structured by host plant 11,15,16,20,25,27,32,34 and heritable bacterial symbionts are involved in expanding diet breadth 26,28 .
Our results show aphids specialized on a single plant family (Asclepiadaceae) that vary in diet breadth have taxonomically similar bacterial communities at the species level (i.e. many taxa in common), but differ in strain diversity and relative abundance of key symbionts (e.g. Arsenophonus). Horizontal transfer of facultative symbionts via host plants can occur in aphids 60,61 , which could contribute to similarities in symbiont communities (i.e. shared ASVs) across different species that overlap in host range and/or naturally co-occur on plants. Another possible explanation is exposure to similar nutritional and chemical profiles 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. Common milkweed (A. syriaca) is among the most abundant species in the midwestern United States 43 . Although A.nerii and A. asclepiadis can colonize multiple milkweed species, it is likely the Indiana populations in this study are primarily exposed to A. syriaca, which could further contribute to similarities in microbial communities. www.nature.com/scientificreports/ Additionally, phylogenetic relatedness can generally result in closely related aphid species harboring more similar microbiota than distantly related species 21 .
Although milkweed aphid microbiomes were overall similar at the species level, bacterial communities may vary across populations and differences in facultative 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 shape-shifting insect symbiont, known best for reproductive manipulation of its host 62,63 . 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 64 . It is also possible Arsenophonus provides a general fitness boost, similar to what has been observed in the soybean aphid 65 . 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 influence each host's symbiotic microbiome 66 , but that volatile organic compounds produced by aphid-associated microbes play a role in attracting ant mutualists 67 . Among milkweed aphids, A.

A. nerii
A. asclepiadis M. asclepiadis    www.nature.com/scientificreports/ asclepiadis is consistently tended by ants and benefits from enhanced protection from predators, while A. nerii is occasionally ant tended and M. asclepiadis appears to be a loner lacking ant friends 68 . 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 influence the intensity of ant attendance in cowpea aphids 69 , it is possible this symbiont has evolved a different ecological role in the case of A. asclepiadis. However, given the 16s rRNA metabarcoding approach used in this study only provides relative abundances, Arsenophonus titer levels will need to be confirmed using additional methods such as quantitative PCR in order to take the first step towards addressing potential functions, including nutritional supplementation or aphid-ant interactions. Detection of additional facultative symbionts (e.g. Serratia) also warrants further investigation into symbiotic relationships and functional roles in milkweed aphid biology and ecology. One unexpected result was the occurrence of a single strain of Buchnera in all three aphid species (i.e. ASV 1; see Fig. 2) even after stringent filtering for false positives. This is the dominant strain infecting A. nerii and is found in significantly higher abundance compared to the other two species (Fig. 2, Supplementary Table S6). The current wealth of research on Buchnera shows this primary endosymbiont of aphids lives intracellularly, relies purely on vertical transmission, and has exhibited co-cladogenesis with aphid hosts over millions of years. It is unlikely this strain has been transferred across species, but it is also unclear why our dataset shows higher levels of sequencing "cross-talk" between samples than previously observed (e.g. 50 ). One possible explanation could be that taxonomic classification using the 16s rRNA gene is unable to provide strain-level resolution for Buchnera in some cases (i.e. multiple strains with identical 16s sequences grouped as ASV1) and therefore additional genomic information is needed to distinguish unique strains found across milkweed aphid species.
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 field. We sampled aphids from a single common host plant rather than multiple milkweed species in order to focus on identifying differences in core heritable symbionts (e.g. presence/absence of taxa, large shifts in relative abundance) and reduce variation introduced by environment and host plant differences. Consequently, our sampling design has limited ability to detect changes in environmentally acquired microbes and does not test for changes in microbiome composition induced by feeding on different host plant species. To gain a more complete understanding of milkweed aphid microbial communities, further work is needed that characterizes microbiota from the broad and narrow specialists (A. nerii & M. asclepiadis) when feeding on multiple host plant species. Imbalance in sampling across aphid species (e.g. fewer M. asclepiadis samples, Supplementary Table S1) resulting from natural variation in prevalence could also mean that some microbial variation was missed. Finally, while this study profiled only bacterial symbionts, additional microbes present in the broader aphid microbiome (e.g. fungi) may be affected by differences in host plant range. Recent work shows bacterial communities associated with milkweed leaves and roots are unique to species, which ultimately shapes the gut microbiomes of insect herbivores like monarch caterpillars that feed on them 70 . Additional studies are therefore needed to dive deeper into the role host plant species plays in shaping milkweed aphid symbiont community composition and function, especially potential links between plant defensive chemistry and microbiome assembly. In general, further research investigating the generalist-specialist gradient using herbivores that feed across multiple plant species and families is needed to clarify the extent to which diet breadth shapes microbial communities (e.g. 24 ).
In summary, we did not find evidence for a gradient in bacterial community diversity associated with variation in diet breadth for milkweed specialized aphid species. Instead, our results suggest overlapping host plant range and shared hosts can result in selection for common microbes and thus microbiomes with high similarity in overall composition when considering the entire pool of microbes detected at the species level. However, milkweed aphids do harbor bacterial populations that vary in strain diversity and relative abundance of Arsenophonus, although a handful of other well-known aphid symbionts were also detected in low abundance. These findings suggest that while diet breadth 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.

Data availability
DNA sequencing data generated and analyzed in this study is available on the NCBI Sequence Read Archive (SRA) repository under accession number PRJNA635683.