Diversity of bacteria associated with Hormaphidinae aphids (Hemiptera: Aphididae).

Bacteria are ubiquitous inhabitants of animals. Hormaphidinae is a particular aphid group exhibiting very diverse life history traits. However, the microbiota in this group is poorly known. In the present study, using high-throughput sequencing of bacterial 16S ribosomal RNA gene amplicons, we surveyed the bacterial flora in hormaphidine aphids and explored whether the aphid tribe, host plant and geographical distribution are associated with the distribution of secondary symbionts. The most dominant bacteria detected in hormaphidine species are heritable symbionts. As expected, the primary endosymbiont Buchnera aphidicola is the most abundant symbiont across all species and has cospeciated with its host aphids. Six secondary symbionts were detected in Hormaphidinae. Arsenophonus is widespread in Hormaphidinae species, suggesting the possibility of ancient acquisition of this symbiont. Ordination analyses and statistical tests show that the symbiont composition does not seem to relate to any of the aphid tribes, host plants or geographical distributions, which indicate that horizontal transfers might occur for these symbionts in Hormaphidinae. Correlation analysis exhibits negative interference between Buchnera and coexisting secondary symbionts, while the interactions between different secondary symbionts are complicated. These findings display a comprehensive picture of the microbiota in Hormaphidinae and may be helpful in understanding the symbiont diversity within a group of aphids.

Introduction

Insects frequently harbor a variety of symbiotic bacteria, and the interactions between them may have important effects on their evolution. Aphids are a group of insects that feed on plant phloem sap and have established a mutualist relationship with the bacterial symbiont Buchnera aphidicola, which inhabits specialized bacteriocytes and supplies essential nutrients that are lacking in the aphid diet (Buchner, 1965; Douglas, 1993; Sandström & Moran, 1999). Buchnera experiences strictly vertical transmission and diversifies parallel to their host during long‐term evolution (Buchner, 1965; Munson et al., 1991; Moran et al., 1993; Baumann et al., 1995; Baumann et al., 1997; Clark et al., 2000; Jousselin et al., 2009; Liu et al., 2013, 2014; Xu et al., 2018).

In addition to the obligate symbiont Buchnera, various secondary symbionts inhabit aphids, namely, Arsenophonus, Fukatsuia symbiotica, Hamiltonella defensa, Regiella insecticola, Rickettsia, Rickettsiella viridis, Serratia symbiotica, Spiroplasma and Wolbachia (Oliver et al., 2010, 2014; Zytynska & Weisser, 2016; Guo et al., 2017). These secondary symbionts distribute erratically in aphids and undergo vertical and some horizontal transmission (Chen & Purcell, 1997; Sandström et al., 2001; Russell et al., 2003; Russell & Moran, 2005; Vorburger et al., 2017; Rock et al., 2018). They play important roles in aphid performance in various environments including protection against parasitic wasps (Oliver et al., 2003, 2005; Vorburger et al., 2010; Hansen et al., 2012; Brandt et al., 2017; Frago et al., 2017), resistance to fungal pathogens (Ferrari et al., 2004; Scarborough et al., 2005; Łukasik et al., 2013), modification of body colors (Tsuchida et al., 2010; Nikoh et al., 2018), interactions with host plants (Leonardo & Muiru, 2003; Tsuchida et al., 2004; Wagner et al., 2015), and thermal tolerance (Chen et al., 2000; Montllor et al., 2002). Furthermore, several facultative symbionts have evolved as co‐obligate endosymbionts to supplement Buchnera, such as Erwinia haradaeae, F. symbiotica, H. defensa, S. symbiotica, Sodalis in some Lachninae species (Pérez‐Brocal et al., 2006; Lamelas et al., 2011; Manzano‐Marín & Latorre, 2014; Manzano‐Marín et al., 2016, 2017, 2019; Meseguer et al., 2017) and Wolbachia in Pentalonia nigronervosa Coquerel (De Clerck et al., 2015; but see Manzano‐Marín, 2019).

The associations between microbial symbionts and aphids varies in different aphid groups. Zytynska and Weisser (2016) reviewed studies about aphid–symbiont associations in 156 aphid species with a strong focus on Western Palearctic samples and mainly Aphidinae and Lachninae groups, 89 and 46 species each, respectively. They revealed that the biological roles of secondary symbionts were dependent on many factors (e.g., aphid species, host plant, genotype); distribution patterns of different symbionts were variable within aphids, which might be contributed by aphid species, host plant species, geography and several environmental factors; and interactions between symbionts were complicated, which might be influenced by both internal (e.g., aphid and symbiont variation) and external factors (e.g., host plant species/abundance, parasitism rate and temperature). Through high‐throughput sequencing, new symbiotic associations with Erwinia‐ and Sodalis‐related bacteria and Type‐X (later named F. symbiotica) were detected in Cinara species (Jousselin et al., 2016; Meseguer et al., 2017). Using the same method, Fakhour et al. (2018) revealed that the compositions of aphid bacterial flora were not limited to commonly known symbionts, and several other bacteria were also present.

However, few studies have revealed the bacterial community and factors that shape it within a large aphid group. Hormaphidinae is an extraordinary group with complex life cycles. Many species in this subfamily are heteroecious, seasonally obligate alternating between primary and secondary plants. They exhibit strong primary host plant specificity, with each tribe feeding on one generic plant, while the associations with secondary host plants are more relaxed, that is, Cerataphidini on Gramineae, Compositae, and Loranthaceae; Hormaphidini on Betula (Betulaceae) and Picea (Pinaceae); and Nipponaphidini on Fagaceae, Lauraceae and Moraceae (Aoki & Kurosu, 2010; Chen et al., 2014). Hormaphidinae aphid species form morphologically diverse galls on primary hosts, secrete a visible wax coating and produce specialized sterile soldiers (Aoki et al., 1977; Aoki & Miyazaki, 1978; Ghosh, 1985, 1988; Stern & Foster, 1996; Chen & Qiao, 2009; Aoki & Kurosu, 2010; Chen et al., 2014). They are mainly distributed in eastern and southeastern Asia (Heie, 1980; Ghosh, 1985, 1988; von Dohlen et al., 2002).

Almost all aphids harbor B. aphidicola as the primary endosymbiont, but in some hormaphidine species, Buchnera has been lost and replaced by yeast‐like symbionts (Fukatsu & Ishikawa, 1992; Fukatsu et al., 1994; Xu et al., 2018). Secondary symbionts in Hormaphidinae have been little studied. Two Hormaphidini species were included in Russell et al. (2003) to explore the distributions of H. defensa, R. insecticola and S. symbiotica, but none of the three symbionts were found in hormaphidine species. Wang et al. (2014) sampled six Hormaphidinae species and detected Wolbachia in all of them. Beyond these studies, the bacterial flora in this extraordinary group remains largely unknown to date.

In the present study, based on extensive taxon sampling, we characterize the microbial communities of Hormaphidinae aphids, evaluate the impact of aphid phylogeny, host plant and geographical distribution on the bacterial community and discuss the symbiont infection patterns and interactions using high‐throughput sequencing of 16S ribosomal RNA (16S rRNA).

Materials and methods

Sampling and extraction of total DNA

Forty‐nine samples representing 23 genera and 49 Hormaphidinae species were collected in this study (Table S1). All specimens were preserved in 95% or 100% and 75% ethanol for molecular experiments and voucher specimens, respectively. All aphid voucher specimens and samples were deposited in the National Zoological Museum of China, Institute of Zoology, Chinese Academy of Sciences, Beijing, China.

Each sample analyzed contained three to 10 individuals from the same aphid clone. Aphid specimens were first immersed in 70% ethanol, washed for 5 min (with vortexing and centrifugation) and then rinsed with sterile water four times to remove body surface contaminations. Total DNA was extracted from pure aphids using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions, and two negative controls were set. The standard cytochrome oxidase subunit I (COI) barcode of each sample was amplified using the primer pair LCO1490/HCO2198 (Folmer et al., 1994) to verify the aphid species identification and to eliminate parasitized aphids. Three other aphid gene sequences, Cytb, EF‐1α and LWO, which were used to reconstruct the phylogeny of Hormaphidinae, were amplified or downloaded from GenBank. All the new sequenced data have been submitted to the GenBank database (Table S1). Phylogenetic congruence between Hormaphidinae species and Buchnera was tested. The more detailed methods of these analyses are provided in the Supporting information, Extended methods.

16S rRNA amplicon amplification and sequencing

The 16S rRNA amplicon of the V3–V4 regions was amplified using the primer pair (341F, 5ʹ‐CCTAYGGGRBGCASCAG; and 806R, 5ʹ‐GGACTACNNGGGTATCTAAT) with a barcode. All polymerase chain reaction (PCR) amplifications were performed in a 30 µL reaction mixture containing 15 µL Phusion® High‐Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, USA), 0.2 µmmol/L forward and reverse primers and approximately 10 ng template DNA. The PCR conditions were as follows: initial denaturation at 98 °C for 1 min, followed by 30 cycles of denaturation at 98 °C for 10 s, annealing at 50 °C for 30 s and elongation at 72 °C for 30 s and final extension at 72 °C for 5 min. Each sample was amplified in duplicate, and one negative and one positive control containing equal amounts of sterile water and DNA of Escherichia coli instead of the aphid DNA were prepared for PCR amplification. The PCR assays for the negative controls (two negative controls in the DNA extraction process and one negative control in the 16S rRNA amplification process), Cerataphis brasiliensis (Hempel), Glyphinaphis bambusae van der Goot and four species in the genus Tuberaphis Takahashi were negative; therefore, these samples were not used for library construction. PCR products were mixed in the same volume with 1× loading buffer (containing SYBR green) and subjected to electrophoresis on a 2% agarose gel. Samples with a bright band between 400 and 450 bp were chosen for further experiments. The target bands of the PCR products were excised and mixed in equidensity ratios and then purified with a GeneJET Gel Extraction Kit (Thermo Fisher Scientific, Massachusetts, United States). A sequencing library was constructed using the NEB Next® UltraTM DNA Library Prep Kit for Illumina (New England Biolabs, Waltham, MA, USA) following the manufacturer's recommendations, and index codes were added. The library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Fisher Scientific) and Agilent Bioanalyzer 2100 system. The library was then sequenced on an Illumina HiSeq2500 platform, and 250 bp paired‐end reads were generated. The raw reads were deposited in the National Center for Biotechnology Information Sequence Read Archive (SRA) database under BioProject accession number PRJNA553318.

Sequencing data analysis

Paired‐end reads were assigned to samples based on their unique barcodes. The reads were then merged using FLASH (V1.2.11) (Magoc & Salzberg, 2011), and low‐quality tags and chimeras were filtered by QIIME (Caporaso et al., 2010). Sequences were clustered into operational taxonomic units (OTUs) at 97% identity by pick_de_novo_otus.py in QIIME. The most abundant sequence was picked as a representative sequence for each OTU to annotate taxonomic information with the RDP Classifier based on the SILVA 132 database (Wang et al., 2007; Quast et al., 2013; Yilmaz et al., 2014). The OTUs of singletons and chloroplasts were excluded (Navas‐Molina et al., 2013). The OTU abundance of each species was rarefied to the value corresponding to the minimum sum of OTU sequences across all the samples to mitigate the differences in the sequencing effort, and then the relative abundance was calculated based on this rarefied abundance data by dividing the abundance of each OTU by the total abundance of a given species. Subsequent diversity analyses were all performed based on this rarefied abundance or relative abundance data. All OTUs assigned to the reported secondary symbionts of aphids were screened out, and the relative abundance of each secondary symbiont was calculated to better explore the symbiont diversity.

Diversity analysis

To evaluate the alpha diversity of aphid bacterial community, the observed species, the Shannon index and the Simpson index of each species were calculated using the phyloseq package in R 3.5.1 (McMurdie & Holmes, 2013; R Core Team, 2018) based on the OTU abundance table. A rarefaction curve was generated based on the index of observed species.

All samples of Hormaphidinae were grouped according to tribe (including three groups) and host plant family (including four groups with samples ≥3) (Table S2). Significance tests of the alpha diversity indices (the Shannon and the Simpson indices) for aphid bacteria and secondary symbionts from different groups were performed using the Kruskal‐Wallis and Tukey's honestly significant difference (HSD) tests implemented in the vegan (Oksanen et al., 2018) and agricolae packages (de Mendiburu, 2017), respectively.

The dissimilarities of the bacterial communities and secondary symbiont communities between samples were quantified by calculating the Bray‐Curtis dissimilarity using the vegan package (Oksanen et al., 2018). The bacterial communities and secondary symbiont communities among groups were clustered using constrained principal coordinate analysis (CPCoA) and nonmetric multidimensional scaling (NMDS) in the vegan package (Oksanen et al., 2018) based on the relative abundance of each genus and the Bray‐Curtis dissimilarity and plots were created in the ggplot2 package (Wickham, 2016). Based on the Bray‐Curtis dissimilarity, permutational multivariate analysis of variance (PERMANOVA or ADONIS) was performed in the vegan package (Oksanen et al., 2018) to discern statistically significant differences as a result of grouping factors, and analysis of similarities (ANOSIM) was used to test whether the dissimilarity between groups was significantly greater than those within groups in the pegas package (Paradis, 2010). Bipartite networks between secondary symbionts and their aphid hosts were constructed based on the relative abundance data using the bipartite package (Dormann et al., 2008). The specificity coefficient (d’) for each secondary symbiont was estimated using the function specieslevel in the bipartite package, which compares the relative abundance of interactions of a secondary symbiont with an aphid species with the average relative abundance of interactions of that particular secondary symbiont across all aphid species (Dormann et al., 2008; Dormann, 2011).

To explore the effect of geographic distance among species on structuring the bacterial community, the Spearman correlation coefficient between beta diversity index (Bray‐Curtis) and geographic distance matrix was tested. A geographic distance matrix was constructed from geographic points (latitudes and longitudes; Table S2) using the GeoDistanceInMetresMatrix function written by Peter Rosenmai. The Spearman correlation coefficient (ρ) between the two matrices was calculated, and the significance of the statistic was evaluated by a permutation procedure using the Mantel test in the vegan package (Oksanen et al., 2018).

Spearman's rank correlation coefficient (ρ) was calculated to explore the interactions between different symbionts associated with Hormaphidinae based on their relative abundance in the Hmisc package (Harrell & with contributions from Charles Dupont and many others, 2018).

Results

Sequencing data

The sequencing of the 16S rRNA V3–V4 amplicons yielded 9 249 188 raw reads. After quality filtering and removal of chimeric sequences, a total of 8 001 564 effective tags with an average length of 427 nt were obtained. The sequences were classified into 3420 OTUs at 97% sequence identity. The rarefaction curve for each sample tended to saturate (data not shown).

Bacterial diversity across Hormaphidinae aphids

After discarding singletons and chloroplast sequences, 3093 OTUs were obtained and annotated to 23 phyla, 203 families and 469 genera. Overall, 43.62% of these OTUs were attributed to Proteobacteria, 18.76% to Firmicutes and 12.58% to Bacteroidetes. The alpha diversity of bacteria in Hormaphidinae was relatively low (mean Shannon index = 0.53, mean Simpson index = 0.25). The bacterial communities were dominated by B. aphidicola, Serratia, Arsenophonus and Wolbachia (Fig. 1). In addition, Gilliamella (the family Orbaceae) was detected in 29 Hormaphidinae species with high relative abundance (average relative abundance across all samples: 1.32%; Fig. S1). The total relative abundance of the above bacteria was more than 93.00% in most samples, and the other bacterial genera accounted for less than 0.50%.

Fig. 1

Bar plot of the top 10 genera within Hormaphidinae. The others includes all relatively low abundance genera excluding the 10 dominant genera.

The primary endosymbiont B. aphidicola was detected in all species with an average relative abundance of 81.73%. Buchnera and the corresponding hormaphidine species are cospeciated. The result of the Jane analysis showed significant phylogenetic concordance between hormaphidine aphids and Buchnera (P < 0.01) (Fig. S2). The ParaFit analysis rejected the null hypothesis that the phylogenetic trees of aphids and Buchnera were randomly associated (ParaFitGlobal = 1.1776, P = 0.0001), and 40 individual host‐parasite‐associated links contributed to the global trace statistic (P < 0.05).

A total of six aphid secondary symbionts were detected in hormaphidine aphids (Fig. 2), and a bipartite network analysis of secondary symbiont interactions with Hormaphidinae species is reported in Figure 3. Arsenophonus inhabited all the species (detection frequency: 43/43; average relative abundance across all samples: 3.79%), followed by Wolbachia (40/43; 2.47%) and S. symbiotica (35/43; 4.15%). H. defensa (11/43; 0.11%) and Rickettsia (21/43; 0.15%) detected in hormaphidine species with low abundance, while R. insecticola was only present in two species with extremely low abundance (<0.004%). Every secondary symbiont contained just several OTUs. There were seven OTUs for Arsenophonus, six for S. symbiotica, four for Wolbachia, two for Rickettsia and only one for H. defensa and R. insecticola. The reads belonging to each OTU were not equal, but there were no more than two OTUs dominating in each aphid species. Noticeably, the same phylotypes were present in some distant aphid species. Almost all of the sampled species (42/43) were infected with at least two secondary symbionts except for Hormaphis betulae Osten‐Sacken, which was only infected with Arsenophonus (Table S3 and S4). The combination of Arsenophonus, S. symbiotica and Wolbachia (12/43) was the most common type, followed by those of Arsenophonus, Rickettsia, S. symbiotica and Wolbachia (10/43) and Arsenophonus, H. defensa, Rickettsia, S. symbiotica and Wolbachia (10/43). Furthermore, Hybothoracaphis laevigata Chen, Jiang, Chen & Qiao was infected with all six detected secondary symbionts. The specificity (d’) of secondary symbionts was inferred (Fig. 3). The specificities of H. defensa, Rickettsia and R. insecticola were 0.57, 0.53 and 0.07, respectively. Arsenophonus, S. symbiotica and Wolbachia had a specificity of zero.

Fig. 2

Phylogeny of hormaphidine aphid species with the relative abundance of the secondary symbionts displayed as bubbles at the tips of the phylogeny. Different colors of the bubbles represent different secondary symbionts and bubble sizes correspond to the relative abundance of each symbiont, as show in the legend.

Fig. 3

Interaction network structure between aphids and secondary symbionts. The width of the links is proportional to the relative abundance of secondary symbionts associated with a given aphid species. The bottom and top boxes represent the Hormaphidinae species and secondary symbionts, respectively. Colors correspond to different secondary symbionts, as shown in the legend. Specificity values (d’) for secondary symbionts are reported.

Comparison of bacterial and secondary symbiont communities associated with Hormaphidinae among different grouping sets

Measurement of within‐sample diversity (alpha diversity) of bacteria showed significant differences between Hormaphidini and both Cerataphidini and Nipponaphidini (Kruskal‐Wallis test: P < 0.05 between Hormaphidini and Cerataphidini and between Hormaphidini and Nipponaphidini for the Shannon and the Simpson index; Tukey's HSD test: P < 0.05 between Hormaphidini and Cerataphidini for the Shannon and the Simpson index; Fig. S3A and B). The microbiota of Hormaphidini including four sampled species had lower diversity than those of Cerataphidini and Nipponaphidini (Fig. S3A and B). However, there were no significant differences in the diversity of secondary symbionts among the three tribes (Fig. 4A and B). The alpha diversity of neither bacteria (Fig. S3C and D) nor secondary symbionts (Fig. 4C and D) showed significant differences among hormaphidine species exploiting different plant families. We found that compositions of the bacterial microbiota and secondary symbiont community were similar among different groups. CPCoA and NMDS did not form any clusters of either the bacterial community (Fig. S3E–H) or the secondary symbiont community (Fig. 4E–H) for three hormaphidine tribes and for species feeding on different plants. Statistical analyses showed that there were no significant differences in bacterial community compositions among the three hormaphidine tribes (ADONIS: F 2,40 = 0.77, R 2 = 0.04, P = 0.57; ANOSIM: R = −0.11, P = 0.97) and among aphids feeding on different plants (ADONIS: F 3,27 = 0.96, R 2 = 0.10, P = 0.42; ANOSIM: R = −0.06, P = 0.68). The composition of the secondary symbiont community was also not different among the three hormaphidine tribes (ADONIS: F 2,40 = 1.17, R 2 = 0.06, P = 0.29; ANOSIM: R = −0.004, P = 0.50) and among aphids feeding on different plants (ADONIS: F 3,27 = 0.83, R 2 = 0.08, P = 0.54; ANOSIM: R = 0.05, P = 0.26). The correlation between the beta diversity index of both bacteria and secondary symbionts and specimen geographic distance was not significant (Mantel test: ρ = −0.08–0.01, P = 0.79–0.39).

Fig. 4

Boxplot comparing the Shannon index of hormaphidine secondary symbiont communities from three tribes (A) and four plant families (C); comparing the Simpson index of hormaphidine secondary symbiont communities from three tribes (B) and four plant families (D). No significant differences across groups on the basis of Kruskal‐Wallis and Tukey's honestly significant difference test. The bottom and top edges of the boxes mark the 25th and 75th percentiles (i.e., first and third quartiles), respectively. Lines within boxes represent the medians, hinges represent the ±25% quartiles, and whiskers represent up to 1.5× the interquartile range. Constrained principal coordinate analysis (CPCoA) plot illustrating the separation of samples based on differences in secondary symbiont community structure among three tribes (5.4% of the total variance, P = 0.29) (E) and four plant families (9.36% of the total variance, P = 0.48) (G). Nonmetric multidimensional scaling (NMDS) plot illustrating the separation of samples based on differences in bacterial community structure among three tribes (stress = 0.07) (F) and four plant families (stress = 0.07) (H). Colors correspond to different groups, as shown in the legend.

Correlation test between different symbionts associated with Hormaphidinae aphids

The results of Spearman's correlation coefficient are shown in Figure 5 and Table S5. The relative abundance of Arsenophonus and Rickettsia (r = 0.48, P < 0.01), Hamiltonella and Regiella (r = 0.44, P < 0.01), Hamiltonella and Rickettsia (r = 0.43, P < 0.01), Hamiltonella and Serratia (r = 0.62, P < 0.001) and Rickettsia and Wolbachia (r = 0.41, P < 0.01) were positively correlated. B. aphidicola had significantly negative correlations with Serratia (r = −0.56, P < 0.001) and Wolbachia (r = −0.57, P < 0.001).

Fig. 5

Heatmap of Spearman correlation coefficients of symbionts. Positive correlations are indicated as red gradients from 0 to 1.0 and negative correlations are indicated as blue gradients from −1.0 to 0, as shown in the legend.

Discussion

Symbionts inhabiting Hormaphidinae aphids

Our study revealed the bacterial communities of Hormaphidinae. B. aphidicola inhabited all the sampled species with the highest relative abundance. Considering the obligate nutritive role of Buchnera and the long‐term endosymbiotic association between aphids and Buchnera (Douglas & Prosser, 1992; Moran et al., 1993; Shigenobu et al., 2000), the ubiquity of the high abundance of Buchnera in the present study seems to be predictable. Buchnera and aphids have been demonstrated to diversify in parallel in several aphid groups (Buchner, 1965; Munson et al., 1991; Moran et al., 1993; Baumann et al., 1995; Baumann et al., 1997; Clark et al., 2000; Jousselin et al., 2009; Liu et al., 2013, 2014; Xu et al., 2018). In the present study, the codiversification of hormaphidine species and the corresponding Buchnera is also confirmed, which presents an instance of parallel evolution between aphids and Buchnera within a subfamily.

Six secondary symbionts were detected in this study; however, the infection patterns of these symbionts varied in Hormaphidinae. Our data showed that Arsenophonus was present in all the sampled Hormaphidinae species with high relative abundance. The ubiquity of Arsenophonus in Hormaphidinae species suggests that this symbiont acquisition could be ancient in this aphid subfamily and followed by vertical transmission. Jousselin et al. (2013) largely surveyed the diversity of Arsenophonus in aphids and revealed a high incidence of Arsenophonus in the Aphis Linnaeus genus; and Zouari et al. (2018) detected Arsenophonus in all studied Aphidini species (two Aphis and three Hyalopterus Koch species). Our results confirm that Arsenophonus is a major secondary symbiont of aphids and is widespread across aphid taxa. Arsenophonus could increase the growth of the soybean aphid population (Wulff & White, 2015), and Aphis craccivora Koch hosting Arsenophonus promoted specialization in locust host plants (Wagner et al., 2015). Similarly, the high prevalence of Arsenophonus in Hormaphidina could not be random, which suggests that this symbiont may play an important role in these species. However, all of these require further investigation.

Wolbachia is a common symbiont of terrestrial arthropods and can manipulate the reproduction of mutualists (Stouthamer et al., 1999; Zug & Hammerstein, 2012). Several studies failed to detect Wolbachia in aphids (West et al., 1998; Tsuchida et al., 2002; Kittayapong et al., 2003; Nirgianaki et al., 2003; Carletto et al., 2008), while Wang et al. (2014) found a widespread infection of Wolbachia in Chinese aphid populations. Consistent with the results of Wang et al. (2014), all but two species in this study hosted Wolbachia. The true figure of Wolbachia diversity might have been underestimated because of its low titer, its high genetic divergence or inappropriate detection methods (Augustinos et al., 2011). However, the effects of Wolbachia in aphids are still unclear, and further detailed studies are needed to illustrate the exact role of Wolbachia.

For aphids suffering from heat shock or subjected to constant high temperature, their fitness increased while hosting S. symbiotica, possibly by rescuing the primary symbiont Buchnera (Chen et al., 2000; Montllor et al., 2002; Russell & Moran, 2006). Field studies showed a higher prevalence of S. symbiotica in aphids collected in the summer season than in aphids collected 2–4 months earlier at the same site (Montllor et al., 2002). The main distribution areas of Hormaphidinae aphids are eastern and southeastern Asia with warm or hot climates (Heie, 1980; Ghosh, 1985, 1988; von Dohlen et al., 2002). In the present study, all samples were collected in subtropical or tropical zones in China. Chronic exposure to high temperature may reduce the fitness of aphids; hence, most sampled species (86.05%) host S. symbiotica, which may increase their thermal tolerance and suitability under high temperature.

The defensive roles of H. defensa, R. insecticola and Rickettsia have been documented in many studies (Ferrari et al., 2004; Oliver et al., 2005; Scarborough et al., 2005; Vorburger et al., 2010; Hansen et al., 2012; Łukasik et al., 2013; Oliver et al., 2014; Vorburger & Rouchet, 2016). However, in the present study, Rickettsia and H. defensa were detected in no more than half of the samples with low abundance (<1%), and R. insecticola was only screened in two of all sampled species with an extremely low abundance (<0.004%). Most species in Hormaphidinae secrete wax to protect against fungal infection, and the visible wax coating can protect against parasitoids and predators (Smith, 1999; Moss et al., 2006; Pope, 2010; Chen & Qiao, 2012; Su et al., 2016). Alternating between different host plants could also help hormaphidine species hide from parasitoids and predators (Way & Banks, 1968; Eastop, 1998). These intrinsic characteristics of Hormaphidinae can reduce the pressures of parasitoids and predators on these aphid groups. Furthermore, carrying defensive symbionts can entail costs in aphids (Chen et al., 2000; Oliver et al., 2008; Polin et al., 2014; Zytynska et al., 2019). Our results revealed negative correlations between the primary symbiont Buchnera and most secondary symbionts, which indicates that competition for limited resources and inhabitancies within aphids between the primary and secondary symbionts may exist (Zytynska & Weisser, 2016). Therefore, the patchy distributions of these symbionts within Hormaphidinae may be a balance selection (Oliver et al., 2014).

We also found a Gilliamella bacterium in more than half of the sampled species in the present study. The relative abundance of this bacterium was high in several species and second only to the primary symbiont Buchnera. Gilliamella is a gut symbiont of bees and stimulates bees to utilize several toxic sugars; therefore, Gilliamella maintains the health of the bee host (Kwong & Moran, 2013; Kwong et al., 2014; Zheng et al., 2016). We propose that this bacterium might also be a gut symbiont in aphids. However, the actual role of Gilliamella in aphids remains to be experimentally tested.

The influence of geographical distribution, host plant and aphid phylogeny on symbiont communities

Geographical distribution and food plants of aphids have been reported to influence symbiont communities, but the samples involved in these studies were mainly different populations from the same species (Najar‐Rodriguez et al., 2009; Jones et al., 2011; Ferrari et al., 2012; Russell et al., 2013; Brady et al., 2014; Henry et al., 2015; Zhao et al., 2016; Gallo‐Franco et al., 2019; Guo et al., 2019). Our results showed that there were no correlations between the community of aphid secondary symbionts and aphid phylogeny or aphid distributions. The aphid food plant also appears to have no impact on the community profiles of bacteria and symbiotic microbes. Neither ordination analyses nor statistical tests revealed effects of both the host plant and hormaphidine tribes on bacteria or secondary symbiont communities. These findings are consistent with previous studies that reported similar results in which the geographical distributions or host plants did not structure symbiont communities (Fakhour et al., 2018), and the presence of certain secondary symbionts was not affected by aphid phylogeny (Henry et al., 2015). No obvious specificities of secondary symbionts toward aphid species were found in the network analysis. Overall, secondary symbionts did not form any specific clusters but showed a relatively uniform distribution across hormaphidine taxa. Furthermore, every secondary symbiont contained just a few OTUs, and the same phylotype was shared by distantly related aphid taxa. These results indicate that horizontal transfers of secondary symbionts may occur in Hormaphidinae. Horizontal transmission has been reported to repeatedly occur in several secondary symbionts (Sandström et al., 2001; Russell et al., 2003; Jousselin et al., 2013). Bacterial symbionts can perform horizontal transfer during aphid sexual reproduction via aphid host plants and through sequential stabbing in different aphids by parasitoids (Moran and Dunbar, 2006; Gehrer & Vorburger, 2012; Chrostek et al., 2017; Pons et al., 2019). Many hormaphidine species are heteroecious holocyclic (Ghosh, 1985, 1988). The species in Hormaphidinae with sexual generation, seasonal host alternation between primary and secondary host plants and repeated migrations among different secondary host plants could greatly increase the possibility of horizontal transmission of their secondary symbionts.

Interactions between secondary symbionts

Multi‐infections with secondary symbionts occurred commonly in Hormaphidinae. Strong positive correlations of various secondary symbiont combinations were revealed by Spearman correlation analysis in this study. The superinfection of secondary symbionts has been documented in several studies (Ferrari et al., 2012; Russell et al., 2013; Smith et al., 2015; Zytynska et al., 2016; Zhang et al., 2019). Multiple infection may be a result of frequent horizontal transmission of secondary symbionts and alternatively form a horizontal gene pool for recombination or transfer (Moran & Dunbar, 2006; Henry et al., 2013; Russell et al., 2013). Coinfections of Hamiltonella–Serratia and Hamiltonella–Fukatsuia exhibited greater resistance to parasites in Acyrthosiphon pisum (Harris) (Oliver et al., 2006; Guay et al., 2009). However, co‐infection with Hamiltonella and Arsenophonus enhanced the self‐fitness of Aphis gossypii Glover rather than the resistance against parasitoids (Ayoubi et al., 2020). In contrast, coinfecting Hamiltonella negatively affected the beneficial phenotype provided by Rickettsiella (Leclair et al., 2017). Furthermore, the cost of hosting multiple symbionts may additively combine (Oliver et al., 2006; Leclair et al., 2017). McLean et al. (2018) showed a polymorphic figure of multiple infections with different symbiont combinations. These findings suggest that the interactions between secondary symbionts can be synergistic, additive or antagonistic.

Conclusions

In this study, using high‐throughput sequencing of bacterial 16S rRNA gene amplicons, we described the bacterial diversity in the aphid subfamily Hormaphidinae. The primary endosymbiont Buchnera unsurprisingly inhabited all species, in accordance with its obligate mutualist role. Otherwise, we provide a good example of codiversification between aphids and Buchnera at the subfamily level. Arsenophonus was the predominant secondary symbiont in Hormaphidinae species, and its high prevalence might indicate an ancient acquisition of this symbiont. There were no relationships between symbiont diversity and any of the aphid tribes, host plants or geographical distributions. These reveal unspecific clusters of secondary symbionts, which suggest horizontal transmission may occur for these secondary symbionts. Moreover, multiple infections of secondary symbionts were common in Hormaphidinae, but the interactions between them were very complicated. In addition, we first reported the bacterium Gilliamella in Hormaphidinae, and the high abundance of Gilliamella indicated that it may exert biological effects on aphids.

Disclosure

The authors declare that they have no conflict of interests.

Table S1. Voucher information and GenBank accession numbers for all aphid species used in this study.

Table S2. Group information.

Table S3. Occurrence of secondary symbionts in Hormaphidinae.

Table S4. Infection pattern of secondary symbiont across Hormaphidinae.

Table S5. Spearman correlation coefficients of symbionts in Hormaphidinae.

Fig. S1. Phylogeny of Hormaphidinae species with the relative abundance of the bacteria Gilliamella displayed as triangle at the tips of the phylogeny. Triangle size corresponds to the different relative abundance according to the inset legend.

Fig. S2. Cophylogeny of Hormaphidinae and Buchnera from Jane 4.0. Blue and black lines indicate the phylogenies of the Buchnera and aphids, respectively. The reconciled trees were from the ML aphid tree and the ML Buchnera tree. Hollow circles indicate cospeciation events, solid circles indicate duplications, solid circles with arrows indicate host switch events, and dashed lines indicate loss events.

Fig. S3. Boxplot comparing the Shannon index of bacterial community among Hormaphidinae samples from three tribes (A) and four plant families (C); comparing the Simpson index of bacterial community among Hormaphidinae samples from three tribes (B) and four plant families (D). Boxes with the same letter are not significantly different, while those without same letters are significantly different on the basis of Kruskal‐Wallis and Tukey's HSD tests. The bottom and top edges of the boxes mark the 25th and 75th percentiles (that is, first and third quartiles), respectively. Lines within boxes represent the medians, hinges represent the +/− 25% quartiles, and whiskers represent up to 1.5x the interquartile range. CPCoA plot illustrating the separation of samples based on differences in bacterial community structure among three tribes (E) and four plant families (G); NMDS plot illustrating the separation of samples based on differences in bacterial community structure among three tribes (F) and four plant families (H). Colors correspond to different groups, as shown in the legend.

Click here for additional data file.