Insights Into the Species-Specific Microbiota of Greenideinae (Hemiptera: Aphididae) With Evidence of Phylosymbiosis.

Aphids and their symbionts represent an outstanding model for studies of insect-symbiont interactions. The aphid microbiota can be shaped by aphid species, geography and host plants. However, the relative importance of phylogenetic and ecological factors in shaping microbial community structures is not well understood. Using Illumina sequencing of the V3-V4 hypervariable region of the 16S rRNA gene, we characterized the microbial compositions of 215 aphid colonies representing 53 species of the aphid subfamily Greenideinae from different regions and plants in China, Nepal, and Vietnam. The primary endosymbiont Buchnera aphidicola and secondary symbiont Serratia symbiotica dominated the microbiota of Greenideinae. We simultaneously explored the relative contribution of host identity (i.e., aphid genus and aphid species), geography and host plant to the structures of bacterial, symbiont and secondary symbiont communities. Ordination analyses and statistical tests highlighted the strongest impact of aphid species on the microbial flora in Greenideinae. Furthermore, we found a phylosymbiosis pattern in natural Greenideinae populations, in which the aphid phylogeny was positively correlated with microbial community dissimilarities. These findings will advance our knowledge of host-associated microbiota assembly across both host phylogenetic and ecological contexts.

Introduction

Sap-feeding aphids harbor a diverse community of microbes, predominantly symbionts (Jousselin et al., 2016; Guyomar et al., 2018), which represent an excellent model to study symbiosis in insects. Almost all aphids possess the primary endosymbiont Buchnera aphidicola within specialized bacteriocytes to compensate for the nutritional deficiency in their diet (Baumann et al., 1995). Buchnera is inherited by strict vertical transmission (Koga et al., 2012) and cospeciates with aphids (Munson et al., 1991; Clark et al., 2000; Xu et al., 2018). In addition, aphids host several types of secondary symbionts that are not essential for survival. The beneficial roles of secondary symbionts have been demonstrated in numerous studies, including protection against parasitoid wasps and fungi (Oliver et al., 2003; Scarborough et al., 2005; Heyworth and Ferrari, 2015), tolerance to heat shock (Chen et al., 2000; Montllor et al., 2002; Guay et al., 2009) and changes in host plant range (Tsuchida et al., 2004; Wagner et al., 2015). Secondary symbionts usually reside in secondary bacteriocytes, sheath cells, or the hemocoel (Fukatsu et al., 2000) and undergo maternal transmission and occasionally horizontal transmission (Russell et al., 2003; Michalik et al., 2014; Pons et al., 2019). To date, nine secondary symbionts have been most frequently studied in aphids, namely, Arsenophonus (Russell et al., 2003), Fukatsuia symbiotica (Guay et al., 2009), Hamiltonella defensa (Darby et al., 2001), Regiella insecticola (Sandström et al., 2001), Rickettsia (Chen et al., 1996), Rickettsiella viridis (Tsuchida et al., 2010), Serratia symbiotica (Unterman et al., 1989), Spiroplasma (Fukatsu et al., 2001), and Wolbachia (Augustinos et al., 2011).

Numerous studies have investigated the ecological factors shaping aphid-symbiotic community structures and highlighted the importance of geography (Tsuchida et al., 2002; Sepúlveda et al., 2017; Gallo-Franco et al., 2019) and host plants (Simon et al., 2003; Ferrari et al., 2012; Brady and White, 2013; Xu S. F. et al., 2020). Symbiont infection patterns were also related to aphid species (Qin et al., 2021a), characteristics of aphids (Xu et al., 2021), and seasonal shifts (Smith et al., 2015; Liu et al., 2019). Considering the impact of host phylogeny, some studies substantiated the correlation between microbiota dissimilarities and aphid relatedness (McLean et al., 2019; Qin et al., 2021a,b), whereas the phylogenetic signature was not detected in some aphid groups (e.g., Eriosomatinae and Hormaphidinae) (Xu T. T. et al., 2020; Xu et al., 2021).

Phylosymbiosis is a term used to describe the eco-evolutionary pattern, in which the microbial community compositional similarity is significantly correlated with host phylogeny (Brucker and Bordenstein, 2013; Lim and Bordenstein, 2020). This pattern does not assume stable and long-term associations between hosts and microbes. Both host-microbe codiversification and ecological filtering by conserved host traits can contribute to the mechanisms underpinning phylosymbiosis (Sanders et al., 2014; Mazel et al., 2018; Lim and Bordenstein, 2020). Evidence for phylosymbiosis has been provided in insects, birds and mammals (Brooks et al., 2016; Groussin et al., 2017; Chiarello et al., 2018; Hammer et al., 2020; Trevelline et al., 2020). Phylosymbiosis signals are variable across host clades and may be weakened or even erased over extended evolutionary periods (Mazel et al., 2018). However, the extent to which host phylogenetic scales (e.g., from intraspecific to family levels) play a role in phylosymbiosis detection remains to be explored.

Greenideinae is a subfamily of Aphididae within the Hemiptera order, in which 179 extant species and 10 fossil species belonging to 18 genera have been recognized worldwide (Favret, 2021). Greenideinae comprises three tribes, namely, Cervaphidini, Greenideini, and Schoutedeniini. Cervaphidini contains one fossil genus and six extant genera including Anomalosiphum and Cervaphis. Greenideini is represented by Allotrichosiphum, Eutrichosiphum, Greenidea, Greenideoida, Mesotrichosiphum, Mollitrichosiphum, and Tritrichosiphum. Schoutedeniini consists of Eonaphis, Paulianaphis, Schoutedenia, and Palaeogreenidea (fossil genus). Greenideinae is monoecious, with a holocyclic or anholocyclic life cycle. Host plants of Schoutedeniini species belong primarily to Euphorbiaceae (Blackman and Eastop, 2021). Cervaphidini colonizes plants from Fagaceae, Fabaceae, Myrtaceae, Tiliaceae, and so on. The majority of Greenideini species feed on young leaves or shoots of Fagaceae. Other plants have also been reported to serve as hosts of Greenideini, such as Betulaceae, Juglandaceae, and Moraceae. Some Greenideini species exhibit high host plant diversity and colonize plants from different families, whereas many species feed on few types of plants. Greenideinae is mainly distributed in eastern and southern Asia, with 92 species recorded in China (Blackman and Eastop, 2021). Some Greenideinae species, such as Greenidea ficicola and Greenidea psidii, are insect pests that seriously threaten the economics of agriculture and horticulture. Qin et al. (2021a) uncovered the phylosymbiotic signatures of the microbial community associated with Mollitrichosiphum at the genus level and revealed intraspecific phylosymbiosis within one aphid species, M. tenicroperus (Qin et al., 2021b). However, little is known about the microbiota structures across the aphid subfamily Greenideinae. Greenideinae provides an opportunity to explore whether phylosymbiotic patterns occur at various host taxonomic levels in a specific aphid group.

Based on high-throughput 16S rRNA gene sequencing, we aimed to reveal the microbial community composition associated with Greenideinae. We fully assessed the contribution of host identity (i.e., aphid genus and aphid species), geographic distribution and host plant on bacterial, symbiont (including Buchnera and secondary symbionts) and secondary symbiont communities. To further understand the ecological and evolutionary processes of microbiota assemblages, we explored the presence of a phylosymbiotic pattern in Greenideinae.

Materials and Methods

Aphid Collection and DNA Extraction

We sampled 215 aphid colonies from 29 plant families and 32 geographic regions of China, Nepal, and Vietnam, representing 53 species and 9 genera within Greenideinae (Table 1). Information regarding sample collection is provided in Supplementary Table 1. The samples for slide mounting were preserved in 75% ethanol, and the specimens for molecular experiments were stored in 95% ethanol at −20°C. All the samples were deposited in the National Animal Collection Resource Center, Institute of Zoology, Chinese Academy of Sciences, Beijing, China. Aphid identification was performed using morphological examination and DNA barcoding.

TABLE 1

Sampling information for the Greenideinae samples used in the present study.

Aphid tribe	Aphid genus	Aphid species	Number of samples
Cervaphidini	Anomalosiphum	A. takahashii, A. tiomanensis	4
	Cervaphis	C. quercus, C. rappardi	5
Greenideini	Allotrichosiphum	A. cyclobalanopsidis	2
	Eutrichosiphum	E. alnicoia, E. alnifoliae, E. apicifuscum, E. dubium, E. heterotrichum,	49
		E. khasyanum, E. kumaoni, E. parvulum, E. pasaniae, E. pseudopasaniae,
		E. sinense, E. sp. 1, E. sp. 2, E. tattakanum
	Greenidea	G. anonae, G. ayyari, G. brideliae, G. bucktonis, G. camelliae,	76
		G. castanopsidis G. cayratiae, G. decaspermi, G. ficicola, G. flacourtiae,
		G. kuwanai, G. nigra, G. nipponica, G. prunicola, G. psidii,
		G. querciphaga, G. sp. 1, G. sp. 2, G. sp. 3, G. sp. 4, G. symplocosis
	Greenideoida	G. longirostrum, G. lutea	2
	Mesotrichosiphum	M. pentaiarticulatum	1
	Mollitrichosiphum	M. luchuanum, M. montanum, M. nandii, M. nigrofasciatum, M. nigrum,	65
		M. rhusae, M. tenuicorpus, M. tumorisiphum
Schoutedeniini	Schoutedenia	S. emblica, S. ralumensis	11

DNA extraction was performed using the whole body of a single adult viviparous female per colony. To remove external microbial contaminants, aphid individuals were initially washed with 70% ethanol for 5 min and then rinsed with sterile water five times. Total DNA was extracted using the DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany). The negative control contained an equal amount of sterile water instead of aphid DNA. To identify aphid species and eliminate parasitized samples, DNA extracts were verified using National Center for Biotechnology Information (NCBI) BLAST searches of the barcode sequence targeting the cytochrome c oxidase subunit I (COI) gene. The DNA samples were kept at −20°C.

16S rRNA Gene Amplification, Illumina Sequencing, and Sequence Processing

The V3–V4 hypervariable region of the 16S rRNA gene was amplified using the primers 341F (5′-CCTAYGGGRBGCASCAG-3′) and 806R (5′-GGACTACNNGGGTATCTAAT-3′). A negative control was set up for the polymerase chain reaction (PCR), and all amplifications were conducted in triplicate. The 30-μL PCR was prepared with 15 μL of Phusion High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, United States), 3 μL primers and 10 ng of template DNA. The PCR thermal program was as follows: 1 min at 98°C for initial denaturation; 30 cycles of 10 s at 98°C for denaturing, 30 s at 50°C for annealing and 30 s at 72°C for elongation; and 5 min at 72°C for final extension.

Amplification products were recovered using 2% agarose gel electrophoresis, and positive samples with a bright band of 400–450 bp were purified with the GeneJET Gel Extraction Kit (Thermo Scientific, Wilmington, DE, United States). Sequencing libraries were constructed with the NEBNext Ultra DNA Library Prep Kit (New England Biolabs, Ipswich, MA, United States) and quantified by a Qubit 2.0 fluorometer (Thermo Scientific, Wilmington, DE, United States) and an Agilent Bioanalyzer 2100 system. Finally, the libraries were sequenced on an Illumina HiSeq 2500 PE250 platform (Illumina, San Diego, CA, United States).

Paired-end reads were assembled with a minimum overlap size of 10 bp and an error rate of 10% using FLASH v1.2.7 (Magoč and Salzberg, 2011). Low-quality sequences were filtered by QIIME v1.9.1 (Caporaso et al., 2010), and chimeras were discarded by UCHIME v4.2.40 (Edgar et al., 2011). Then, sequences were clustered into one operational taxonomic unit (OTU) at 97% similarity. Taxonomic annotations of OTUs were assessed with the RDP classifier (Wang et al., 2007) against the SILVA 128 reference database (Quast et al., 2013). We also manually checked the taxonomic affiliation using BLAST searches against the GenBank database. The OTUs with abundances below 0.005% were filtered as described by Bokulich et al. (2013). The mean sequence number of three PCR replicates per sample was calculated to estimate the abundance of each OTU. To mitigate the heterogeneity of sequencing depth, each sample was rarefied to the minimum value across all the samples (53,497 reads) using the ‘otutab_norm’ function in USEARCH v10.0 (Edgar, 2010). Finally, we obtained an OTU table containing taxonomic classifications of bacterial OTUs and the sequence number per sample (Supplementary Table 2a).

Statistical Analyses

To better explore the variations in microbial community structures associated with Greenideinae, two reduced OTU tables, including symbionts (Buchnera and secondary symbionts) (Supplementary Table 2b) and secondary symbionts (Supplementary Table 2c), of aphids were generated. The relative abundances of OTUs were calculated using the decostand function and total method in the R package vegan (Oksanen et al., 2018). The distribution and relative abundances of symbionts were mapped against the phylogeny of Greenideinae aphids (detailed analysis methods are provided in the Supplementary Material) using iTOL v6.3.2 (Letunic and Bork, 2021). We implemented all of the subsequent statistical analyses with the bacterial, symbiont and secondary symbiont OTU tables. Samples were grouped according to aphid genus, aphid species, geographic distribution and host plant. Detailed grouping information is shown in Supplementary Table 3. Statistical analyses were conducted on all groups and groups with a sample size ≥ 3, excluding Mantel tests, Procrustes analyses, and ancestral secondary symbiont reconstruction. Samples with ambiguous information were excluded from analyses.

Shannon and Simpson indices were assessed to evaluate the alpha diversity using the diversity function of vegan. Considering the non-normal distribution (Shapiro–Wilk test, p < 0.05) of alpha diversity data, we used non-parametric Kruskal–Wallis tests to investigate the microbiota variation among groups. To simultaneously estimate the separate contributions of aphid genus (or aphid species), geography and host plant to microbial alpha diversity, three-way analysis of variance (ANOVA) was performed using the avop function in the lmPerm package (Wheeler and Torchiano, 2010). Three-way ANOVA was conducted on groups with sample size ≥ 2 and sample size ≥ 3, as the groups with one sample were insufficient for this analysis.

The Bray–Curtis distance (vegdist function of vegan) and unweighted UniFrac distance (GUniFrac function of GUniFrac) (Chen and Chen, 2018) were calculated to quantify beta diversity. First, we applied two kinds of ordination analyses, namely, unconstrained non-metric multidimensional scaling (NMDS) and constrained principal coordinate analysis (cPCoA), to visualize the dissimilarity of microbial communities based on the beta diversity. Based on Bray–Curtis and unweighted UniFrac distance matrices, NMDS (stress values < 0.05 indicate excellent representation) was performed with the metaMDS function in vegan. CPCoA (capscale and anova.cca functions in vegan) was implemented with Bray–Curtis distance, as unweighted UniFrac distance was not suitable for this analysis. Then, the statistical significance of the microbiota variation between groups was determined by permutational multivariate analysis of variance (PERMANOVA) and analysis of similarities (ANOSIM). PERMANOVA and ANOSIM were conducted on Bray–Curtis and unweighted UniFrac distances using the adonis and anosim functions, respectively, with 10,000 permutations in vegan.

To explore the relationship between aphid phylogeny and microbiota dissimilarities, the Mantel test (mantel function) and Procrustes analysis (procrustes and protest functions) were performed with all samples using the matrices of aphid cophenetic distances and beta diversity (Bray–Curtis and unweighted UniFrac distances) in vegan. Cophenetic distances calculate the pairwise distances between the pairs of tips from a phylogenetic tree using its branch lengths. We used the cophenetic function in the package ape (Paradis and Schliep, 2019) to calculate the cophenetic distances of aphid samples. The Mantel test is a commonly used approach to quantify the correlation between two matrices (Anderson and Walsh, 2013). Procrustes analysis is more powerful, in which two matrices are scaled and rotated to maximize their similarity using principal component analysis (PCA) (Peres-Neto and Jackson, 2001).

To investigate the microbial composition changes that occurred during the evolutionary history of Greenideinae, a presence/absence matrix of six secondary symbionts in this study was mapped onto an aphid phylogenetic tree. The relationships of aphid species were simplified from the maximum-likelihood tree of Greenideinae aphid samples. Ancestral microbiome reconstruction was performed using parsimony and Bayesian approach. Parsimony reconstruction was implemented in Mesquite v3.70 (Maddison and Maddison, 2021) with the ‘trace character history’ option and unordered character state transformations. For Bayesian reconstruction, we applied the reverse jump Markov chain Monte Carlo framework with a prior drawn from a uniform distribution in BayesTraits V3.0.5 software (Pagel and Meade, 2006). Three independent analyses were performed for a total of 6,000,000 generations, sampling every 1000 generations. The first 25% of the generations were removed as burn-in, and the acceptance rates were between 0.2 and 0.4.

Results

Microbial Community Composition Across Greenideinae Aphids

After all quality filtering, a total of 11,506,428 reads (53,518 reads per sample) were obtained. The 16S rRNA gene sequences were clustered into 162 OTUs and annotated to 35 genera, 24 families, 16 orders, 12 classes, and 6 phyla of bacteria. At the phylum level, the microbial community of Greenideinae was composed mainly of Proteobacteria, with an average relative abundance of 98.74%. Gammaproteobacteria (95.15%) and Enterobacteriales (94.72%) were the most commonly classified class and order, respectively. The most abundant family was Enterobacteriaceae (94.42%), followed by Anaplasmataceae (2.10%) and Rickettsiaceae (1.34%). At the genus level, seven symbionts were detected, among which the relative abundances of Buchnera (68.19%), S. symbiotica (14.97%), Wolbachia (2.10%), Arsenophonus (1.80%), and Rickettsia (1.34%) were greater than 1% (Supplementary Table 4).

The primary endosymbiont Buchnera was present in all the samples and predominated in most. In addition, each examined sample was simultaneously infected with 2–6 secondary symbionts. S. symbiotica and Wolbachia were the most common, with an infection frequency of 215/215, followed by Arsenophonus (171/215), Rickettsia (96/215), H. defensa (88/215), and F. symbiotica (86/215). S. symbiotica was the most abundant secondary symbiont in Greenideinae, and its relative abundance was even higher than that of Buchnera in the aphid genus Schoutedenia. The infection pattern of secondary symbionts within Greenideinae was variable among different aphid genera and species (Figure 1). The most frequent secondary symbiont infection type was the combination of S. symbiotica, Wolbachia, and Arsenophonus (50/215) (Supplementary Table 5).

FIGURE 1

Symbiont community composition shown on the maximum-likelihood tree of Greenideinae aphids. Gray dots on the phylogeny nodes indicate bootstrap support of over 50%.

At the OTU level, Buchnera (Supplementary Figure 1) and S. symbiotica (Figure 2) harbored 20–30 OTUs, whereas other secondary symbionts were represented by only 1–6 OTUs (Figure 3 and Supplementary Figures 2, 3). OTU11 or OTU 4915 of Buchnera predominated in most samples of the aphid genera Cervaphis and Mollitrichosiphum (Supplementary Figure 1). The aphid genera Greenidea and Eutrichosiphum were characterized mainly by more than one type of OTU belonging to Buchnera. Regarding the secondary symbionts, the dominant OTUs of S. symbiotica usually differed among aphid species, although most OTUs were widely distributed in Greenideinae (Figure 2). Wolbachia was represented by four OTUs, among which the relative abundance of OTU4 was high in some samples (Figure 3). We also observed the variability of the relative abundance and diversity of OTUs belonging to Arsenophonus associated with Greenideini species (excluding aphid species of Anomalosiphum, Cervaphis, and Schoutedenia) (Supplementary Figure 2). The prevalence and relative abundance of other secondary symbionts (i.e., Rickettsia, H. defensa, and F. symbiotica) were quite low in most samples (Supplementary Figure 3). For example, H. defensa and Rickettsia were not detected in the tribe Cervaphidini, which exhibited the lowest secondary symbiont diversity in the present study. The majority of Greenideinae aphids rarely host Rickettsia, except samples from Mollitrichosiphum.

FIGURE 2

Heatmap representing the distribution and relative abundances of Serratia symbiotica OTUs among Greenideinae aphids. The phylogenetic relationships of Greenideinae and S. symbiotica OTUs based on maximum-likelihood analyses are presented.

FIGURE 3

Heatmap representing the distribution and relative abundances of Wolbachia OTUs among Greenideinae aphids. The maximum-likelihood trees display the phylogenetic relationships of Greenideinae and Wolbachia OTUs.

Variations in the Microbial Community Diversity of Greenideinae

Kruskal–Wallis tests of alpha diversity revealed a significant variation in the bacterial, symbiont and secondary symbiont communities among the aphid genera, aphid species and host plants (P < 0.05 for both Shannon and Simpson indices). A significant effect of geography on bacterial and symbiont communities was detected (Shannon, P = 0.014–0.040; Simpson, n ≥ 3, P = 0.021–0.027), except that there was no significant effect on the Simpson index upon using all samples (P = 0.057–0.066). However, the secondary symbiont communities did not differ significantly among geographic regions (P = 0.276–0.538).

Three-way ANOVA for alpha diversity indices highlighted the contribution of aphid identity (i.e., aphid genus or aphid species) to shaping the microbial community structures when the effect of different factors was estimated simultaneously. We found that the microbial communities differed significantly among aphid genera [Supplementary Table 6; n ≥ 2, F(7,70) = 2.667–4.292, P ≤ 0.013; n ≥ 3, F(4,69) = 5.333–7.833, P < 0.001] using all types of data in three-way ANOVA, whereas the effect of the host plant was not significant (P = 0.127–0.636). A significant effect of geography was found only on bacterial and symbiont communities with sample size ≥ 3 compared with the effect of aphid genus and host plant [F(14,69) = 1.810–2.000, P = 0.005–0.027]. The results of three-way ANOVA also showed that the bacterial and symbiont communities were primarily structured by aphid species [Supplementary Table 7; P < 0.001; n ≥ 2, F(28,38) = 9.333–11.375; n ≥ 3, F(15,33) = 18.625–19.600]. Aphid species, geography and host plant did not have a significant impact on the secondary symbiont community (P = 0.174–0.924).

Unconstrained NMDS plots based on all types of beta diversity matrices did not show meaningful clustering of samples structured by aphid genus, aphid species, geographic region or host plant (Supplementary Figures 4–7). Conversely, cPCoA plots of Bray–Curtis distances uncovered the significant patterns constrained by these four factors in bacterial, symbiont and secondary symbiont communities (Figure 4 and Supplementary Figure 8; P = 0.001). Aphid species (33.4–62.3% of variance) explained more of the overall variance in the data than aphid genus (9.19–32% of variance), geographic region (13–18.1% of variance) and host plant (12.7–26.9% of variance) in all the analyses.

FIGURE 4

Constrained principal coordinate analysis (cPCoA) plots of Bray–Curtis distances of bacterial (A,D,G,J), symbiont (B,E,H,K), and secondary symbiont (C,F,I,L) communities (n ≥ 3). Plots are structured by aphid genus (A–C), aphid species (D–F), geographic region (G–I), and host plant (J–L).

PERMANOVA (R2 = 0.408–0.725, P < 0.001) and ANOSIM (R = 0.278–0.752, P < 0.001) corroborated the greatest effect of aphid species on structuring the microbial communities based on Bray–Curtis and unweighted UniFrac distances with sample size ≥ 3, although no statistical significance among aphid species (PERMANOVA: P = 0.108–0.141; ANOSIM: P = 0.127–0.391) or aphid genera (PERMANOVA: P = 0.726–0.778; ANOSIM: P = 0.386–0.548) was observed in most cases based on Bray–Curtis distances using all samples (Table 2). The significant impacts of geographic region (P ≤ 0.006) and host plant (P ≤ 0.047) were revealed by PERMANOVA using all types of beta diversity data. However, their contributions (geographic region: R2 = 0.145–0.643; host plant: R2 = 0.113–0.319) were usually limited compared to that of aphid species (R2 = 0.262–0.725). ANOSIM suggested significant dissimilarities in microbial communities among geographic regions using all samples (R = 0.325–0.576, P < 0.001). A significant effect of the host plant was found only in analyses of the Bray–Curtis distance (R = 0.135–0.177, P ≤ 0.001).

TABLE 2

Results of ANOSIM and PERMANOVA based on Bray–Curtis and unweighted Unifrac distances.

Beta diversity distance	Microbial community	Sample size	Aphid genus		Aphid species		Geographic region		Host plant
			ANOSIM (R, P)	PERMANOVA (R2, P)	ANOSIM (R, P)	PERMANOVA (R2, P)	ANOSIM (R, P)	PERMANOVA (R2, P)	ANOSIM (R, P)	PERMANOVA (R2, P)
Bray–Curtis	Bacteria	n ≥ 1	0.001, 0.461	0.033, 0.726	0.030, 0.142	0.261, 0.108	0.349, <0.001	0.420, <0.001	0.135, <0.001	0.311, <0.001
		n ≥ 3	0.714, <0.001	0.438, <0.001	0.752, <0.001	0.713, <0.001	0.025, 0.165	0.150, <0.001	0.138, <0.001	0.214, <0.001
	Symbionts	n ≥ 1	0.004, 0.386	0.032, 0.778	0.031, 0.127	0.260, 0.141	0.325, <0.001	0.433, <0.001	0.139, <0.001	0.319, <0.001
		n ≥ 3	0.685, <0.001	0.466, <0.001	0.713, <0.001	0.725, <0.001	0.008, 0.358	0.151, 0.001	0.148, <0.001	0.224, <0.001
	Secondary symbionts	n ≥ 1	–0.003, 0.548	0.035, 0.728	0.009, 0.391	0.262, 0.024	0.374, <0.001	0.334, <0.001	0.145, 0.001	0.253, <0.001
		n ≥ 3	0.297, <0.001	0.158, <0.001	0.555, <0.001	0.487, <0.001	0.033, 0.145	0.145, <0.001	0.177, <0.001	0.156, <0.001
Unweighted Unifrac	Bacteria	n ≥ 1	–0.010, 0.677	0.041, 0.293	–0.013, 0.652	0.256, 0.203	0.576, <0.001	0.519, <0.001	0.042, 0.187	0.246, <0.001
		n ≥ 3	0.377, <0.001	0.270, <0.001	0.469, <0.001	0.442, <0.001	0.085, 0.008	0.183, <0.001	0.058, 0.122	0.146, <0.001
	Symbionts	n ≥ 1	0.057, 0.002	0.052, 0.080	0.047, 0.056	0.292, 0.014	0.520, <0.001	0.643, <0.001	–0.068, 0.966	0.207, 0.006
		n ≥ 3	0.423, <0.001	0.377, <0.001	0.399, <0.001	0.530, <0.001	0.031, 0.139	0.151, 0.005	–0.063, 0.939	0.116, 0.001
	Secondary symbionts	n ≥ 1	0.053, 0.004	0.062, 0.030	0.066, 0.016	0.327, 0.001	0.415, <0.001	0.631, <0.001	–0.072, 0.978	0.208, 0.047
		n ≥ 3	0.318, <0.001	0.281, <0.001	0.278, <0.001	0.408, <0.001	0.013, 0.328	0.156, 0.006	–0.052, 0.904	0.113, 0.004

Statistically significant P-values (P < 0.05) are highlighted in italics.

Relationship Between Microbial Communities and Aphid Relatedness

Using the Mantel test with all types of beta diversity data, we detected a positively significant correlation between aphid phylogeny and the microbial profiles of bacterial, symbiont and secondary symbiont communities (Table 3; r = 0.145−0.555, P < 0.001). The results of Procrustes analysis were consistent, which showed that the microbial community structures were significantly related to aphid phylogeny (Figure 5; M2 = 0.570−0.938, P = 0.001).

TABLE 3

Correlations between microbial beta diversity and aphid phylogeny estimated by the Mantel test.

Microbial community	Beta diversity distance	Mantel test
		r	P
Bacteria	Bray–Curtis	0.555	<0.001
	Unweighted Unifrac	0.234	<0.001
Symbionts	Bray–Curtis	0.531	<0.001
	Unweighted Unifrac	0.249	<0.001
Secondary symbionts	Bray–Curtis	0.243	<0.001
	Unweighted Unifrac	0.145	<0.001

Statistically significant P-values (P < 0.05) are highlighted in italics.

FIGURE 5

Procrustean superimpositions for PCA-scaled aphid cophenetic distances vs. variations in bacterial (A,D), symbiont (B,E), and secondary symbiont (C,F) communities. Bray–Curtis (A–C) and unweighted UniFrac distances (D–F) were used to estimate the microbiota variations.

Ancient Aphid-Secondary Symbiont Associations

The ancestral microbiome state was reconstructed to further understand shifts in the secondary symbionts over the evolutionary history of Greenideinae. Both parsimony and Bayesian analyses supported the acquisition of S. symbiotica and Wolbachia at the base of Greenideinae (Supplementary Figure 9; Figure 6, and Supplementary Table 8). Regarding the associations of other secondary symbionts and Greenideinae aphids, some incongruence was observed between different analysis approaches. The parsimony-based inference suggested that Arsenophonus was acquired by the most recent common ancestor of Greenideinae and subsequently lost in several aphid species. In addition, parsimony reconstruction reflected the separate gains of secondary symbionts in different aphid clades, such as Rickettsia in the clade containing Mesotrichosiphum, Schoutedenia, Greenideoida, Allotrichosiphum, Eutrichosiphum, and Mollitrichosiphum (Supplementary Figure 9). The results of Bayesian ancestral state reconstruction are summarized in Supplementary Table 8 and visualized in Figure 6; the results indicated two independent acquisition events for H. defensa in the clade containing Eutrichosiphum and Mollitrichosiphum.

FIGURE 6

Ancestral associations of secondary symbionts and Greenideinae aphid species estimated by Bayesian reconstruction. Pie charts at nodes and tips show the presence of different secondary symbionts. The mean posterior probabilities at nodes are given in Supplementary Table 8.

Discussion

Symbiont Diversity Associated With Greenideinae Aphids

The microbial community composition of Greenideinae was dominated by aphid symbionts, among which seven symbionts were identified in the ten most abundant genera. The ubiquity and relative abundance of Buchnera substantiated its obligate nutritive role in the long-term endosymbiotic relationship with host aphids (Munson et al., 1991; Douglas, 1998; Liu et al., 2013; Xu et al., 2018). Multiple infections (i.e., infections with more than one type of symbiont in a host individual) of secondary symbionts were common in Greenideinae aphids. Positive interactions related to cohabitation of secondary symbionts have been demonstrated in many studies, such as higher resistance to parasitism than to single infection (Oliver et al., 2006). However, some coinfecting secondary symbionts with no conditional mutualism can persist in aphids by hitchhiking with beneficial symbionts (Smith et al., 2015).

Serratia symbiotica and Wolbachia were found in all examined Greenideinae samples with high relative abundance. Our study confirmed the widespread distribution of Wolbachia in aphids reported by Wang et al. (2014). It is worth noting that the relative abundance of S. symbiotica was higher than that of Buchnera in all samples of the aphid genus Schoutedenia. Rübsaamen (1905) has documented that Schoutedenia can induce galls on its host plants, differing greatly from other aphids in Greenideinae. This distinctive biological trait may give rise to different symbiont community structures and even novel nutrition-biosynthetic pathways provided by symbionts in the host Schoutedenia. The contribution of S. symbiotica to nutritional biosynthesis has been reported in some aphids (Lamelas et al., 2011; Mccutcheon and Moran, 2012; Bennett and Moran, 2015). Further research is required to illustrate the exact effects of S. symbiotica associated with Schoutedenia aphids.

Arsenophonus associated with Greenideini species exhibited various infection patterns at the OTU level. Liu et al. (2015) suggested that the acquisition of new host plants might account for species differentiation in Greenideini. Improvement of performance on novel host plants and other ecological benefits have been demonstrated in previous studies of Arsenophonus (Wagner et al., 2015; Lenhart and White, 2020). However, further experiments should be performed to determine whether Arsenophonus infection is related to the broadening of the host plant range in Greenideini.

Other secondary symbionts (i.e., Rickettsia, H. defensa, and F. symbiotica) were not frequent symbiotic partners, and their relative abundances were quite low in most samples. Aphid species in Greenideinae usually move rapidly and have long siphunculi that enable them to efficiently emit alarm pheromones (Mondor et al., 2002). In addition, ant attending was observed in many Greenideinae species (Blackman and Eastop, 2021), providing defense from natural enemies for aphids. Prior studies have reported the trade-off between fitness costs (e.g., reduced lifespan and fecundity) and benefits of hosting symbionts (Oliver et al., 2014; Leybourne et al., 2020; Zytynska et al., 2021). Therefore, aphids generally tend not to host defensive symbionts when they obtain protection from ecological habits (Henry et al., 2015). Greenideinae may have defensive effects derived from these life history traits, reducing the maintenance of secondary symbionts providing similar resistance to parasitism, such as Rickettsia, H. defensa, and F. symbiotica.

Contribution of Different Factors to Shaping the Microbiota of Greenideinae

The results of alpha and beta diversity analyses highlighted the role of host identity in shaping the structures of bacterial, symbiont, and secondary symbiont communities associated with Greenideinae. We revealed the strongest impact of aphid species on microbial profiles, which was in line with previous evidence for the existence of a species-specific microbiota in Mollitrichosiphum aphids (Qin et al., 2021a). In addition, significant microbial community structures constrained by the aphid genus were identified. Acquiring a large amount of species-level taxonomic information on insect specimens is usually time-consuming and challenging. Our results suggest that genus-level insect identification is meaningful for investigating host-associated microbiota structures with limited taxonomic knowledge.

Prior studies have demonstrated that geography (Tsuchida et al., 2002; Guo et al., 2019; Xu S. F. et al., 2020) and host plant (Brady et al., 2014; Gauthier et al., 2015; Xu T. T. et al., 2020) are important factors structuring the symbiont community composition associated with aphids. Distinct environmental conditions (Sepúlveda et al., 2017) and spatial limitations on bacterial dispersal (Moeller et al., 2017) can lend to the variability of symbionts. In the present study, the geography and host plant significantly influenced the microbial profiles across Greenideinae, although their contributions were usually weaker than that of host identity. Symbiont-assisted speciation has been reviewed by Brucker and Bordenstein (2012). Some secondary symbionts of aphids, such as Arsenophonus (Wagner et al., 2015), play a crucial role in the exploitation and specialization of host plants. We observed OTU variability in Arsenophonus among Greenideini species. The species differentiation of Greenideini was involved in the enhancement of host plant utilization (Liu et al., 2015). Hence, it seems likely that secondary symbionts mediate dietary breadth in Greenideini, which potentially facilitates the diversification of aphids.

Phylosymbiosis Between Greenideinae Aphids and Their Microbiota

The pattern of phylosymbiosis revealed in Greenideinae was characterized by a positively significant correlation between aphid phylogeny and microbial community dissimilarities. Mantel tests and Procrustes analyses using all types of data consistently showed the phylosymbiotic microbiota assemblage in bacterial, symbiont and secondary symbiont communities. Phylosymbiosis has been reported in host taxa with evolutionary histories ranging from approximately 0.3–108 million years (Brooks et al., 2016). However, several factors are likely to erode phylosymbiosis signals in natural populations during long evolutionary periods, such as horizontal transmission of microbes, complete dietary transition and changes in environmental filters. Combined with previous studies of Mollitrichosiphum (Qin et al., 2021a) and M. tenicroperus (Qin et al., 2021b), our studies revealed that phylosymbiosis can be observed at the subfamily, genus and species levels within natural populations of specific aphid groups, spanning recent host speciation events to more distant host divergence.

Considering the strictly maternal inheritance of Buchnera and primarily vertical transmission of some secondary symbionts in aphids, one of the possible mechanisms underlying phylosymbiosis is evolutionary processes, such as codiversification of Greenideinae and symbionts. The ancestral association between Greenideinae and two secondary symbionts, S. symbiotica and Wolbachia, was supported by ancestral microbiome reconstruction based on both parsimony and Bayesian approaches. The diversification tendency of S. symbiotica at the OTU level also indicated host-symbiont codiversification within Greenideinae. In this study, Buchnera, S. symbiotica and Wolbachia may serve as hub microbes, allowing determination of the whole microbial community composition through microbe–microbe interactions (Fisher and Mehta, 2014; Agler et al., 2016). Furthermore, previous studies did not detect a correlation between microbiota similarities and aphid phylogeny in heteroecious Eriosomatinae (Xu T. T. et al., 2020). The frequent horizontal transmission of secondary symbionts might have erased the phylosymbiosis signals in this aphid subfamily. In contrast, the monoecious life cycle of Greenideinae could have reduced the horizontal transmission of secondary symbionts among aphid populations colonizing different host plants. Therefore, codiversification of key symbionts and relatively stable microbial community composition may be responsible for the phylosymbiosis of Greenideinae.

Ecological filtering by conserved host traits (Moran and Sloan, 2015; Mazel et al., 2018) can also contribute to the mechanisms underpinning phylosymbiosis. Closely related hosts tend to harbor a similar microbiota owing to similar physiologies or immune responses regulating microbes. In addition, diet is an important filter shaping the host-associated microbiota, and complete dietary shifts over long evolutionary periods can even disrupt phylogenetic signals in microbial communities (Muegge et al., 2011; Groussin et al., 2017; Chiarello et al., 2018). Beta diversity analyses showed the influence of host plants on the structures of the microbial communities associated with Greenideinae. Dietary shifts, including direct switching to novel host plants and expansion of the host plant range, did not cover the phylosymbiosis signals in the present study. Host plants might be filters that are phylogenetically correlated with aphids, considering their promotion of species differentiation in Greenideinae. It seems reasonable that filtering by host plants is another driver yielding the phylosymbiotic microbiota of Greenideinae.

Conclusion

Unraveling the influence of ecological and evolutionary factors on microbial community assembly is crucial to understanding symbiosis in nature. We provided the first systematic microbiota landscape associated with Greenideinae. The significant impacts of host identity, geography and host plant on microbial community structures were uncovered, among which the primary contribution of aphid species was highlighted. Moreover, our study found that microbial community dissimilarities are correlated with aphid relatedness, providing evidence for phylosymbiosis in natural aphid populations. We propose that phylosymbiosis in Greenideinae relies on multiple mechanisms. Aphid-symbiont codiversification and filtering by host plants might contribute to the phylosymbiotic microbiota in Greenideinae.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GenBank under accession numbers OL619437–OL619530, OL624544–OL624640, and OL631781–OL631898. Raw 16S rRNA gene amplicon reads were deposited in the NCBI Sequence Read Archive under BioProject accession number PRJNA783284 and PRJNA637573.

Author Contributions

GQ and LJ designed the project. GQ, JC, and LJ identified voucher specimens. MQ conducted molecular experiments and all analyses and wrote the manuscript. JC assisted data analyses. All the authors contributed to revisions.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.