Low bacterial community diversity in two introduced aphid pests revealed with 16S rRNA amplicon sequencing

Sample collection and DNA extraction

A total of 80 individuals of the aphid S. avenae and 52 individuals of the aphid R. padi were collected from oat (Avena sativa) and wheat (Triticum aestivum) crops in two different agroclimatic regions (Maule and Los Ríos regions) in Chile (Table 1). In addition, the field experiments performed in this study were approved for Ethical scientific committee of the Universidad de Talca in Chile (FONDECYT project 3140299). DNA extraction was individually performed for each aphid specimen using the “Salting out” method described by Sunnucks & Hales (1996). The quantification and quality of the extracted DNA was examined by absorbance using Infinite 200 PRO NanoQuant (TECAN) and by electrophoresis in 0.8% agarose gels. Each individual DNA extraction was normalized to a concentration of 5 ng/ul and kept at −20 °C until later 16 S library preparation.

16S rRNA amplicon sequencing library preparation

In order to produce DNA pools that represent the genetic diversity of aphids from different species, locations and host-plants; four DNA pools of 20 S. avenae aphids and six DNA pools of 9 R. padi aphids were used (Table 1). Pools of the genomic DNA were generated in two steps using the Illumina MiSeq protocol for 16S amplicon sequencing (https://support.illumina.com/downloads/). The V3 and V4 variable regions of the 16S ribosomal RNA gene were amplified using the primer pair S-D-Bact-0341-b-S-17 and S-D-Bact-0785-a-A-21 (Klindworth et al., 2012); these regions have a total length of approximately 460 bp. Specific Illumina Adapter sequences were added to the 5′ region of the forward (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG) and reverse primers (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG). PCR reactions were performed in a volume of 15 µl; 1.5 µl of normalized DNA (5 ng/µl), each primer at 0.2 µM plus 7.5 µl of 2x KAPA Hifi HotStart Ready Mix (KAPA Biosystems, Wilmington, MA, USA). The PCR program consisted of initial denaturation at 95 °C for 3 min, followed by 25 cycles of:  95 °C for 30 s, 55 °C for 30 s and 72 °C for 30 s, and then a final extension of 72 °C for 5 min. The amplicons were purified using a PCR clean-up protocol including 10 µl of PCR product, 8 µl of AMpureXP (Beckman Coulter, Brea, CA, USA), and 200 µl of 80% ethanol. The mixture was incubated on a magnetic stand and diluted in 22.5 µl of Tris (10 mM), pH 8.5. The expected size of PCR amplicons was verified using the DNF-900 Kit for Fragment Bioanalyzer (Advanced Analytical Technologies, Inc., Ankeny, IA, USA); the quantity of amplicons was estimated by fluorescence using the Quant-iT PicoGreen dsDNA Assay kit (Invitrogen, Carlsbad, CA, USA) and a HOEFER DQ300 fluorometer (Hoefer Inc., Holliston, MA, USA). Each DNA pool was constructed by mixing 5 µl of each of the amplicons (Table 1). Then DNA pools were subjected to a second PCR where dual indices and Illumina sequencing adapters were attached using a NEXTERA XT index Kit (Illumina, San Diego, CA, USA). This second PCR was conducted in a total volume of 50 µl which contained 5 µl of each pooled DNA, 5 µl of each Nextera XT index primer, 25 µl of 2x KAPA Hifi HotStart Ready Mix, and 10 µl of PCR grade water. PCR program consisted of initial denaturation at 95 °C for 3 min, followed by eight cycles of: 95 °C for 30 s, 55 °C for 30 s and 72 °C for 30 s, and 72 °C for 5 min. The PCR product was corroborated using a Fragment Analyzer and the DNF 479 kit. Finally, each DNA pool was normalized to a concentration of 4 nM and then pooled. The mix DNA pool was prepared for sequencing using the Denature and Dilute Libraries Guide. Paired-end sequencing was performed using the Miseq Reagent Kit v3 (2 × 300 cycles) on the MiSeq Illumina sequencing platform in the AUSTRAL-omics Core-Facility (Facultad de Ciencias, Universidad Austral de Chile).

Data analysis

Removal of adapters and quality filtering of the data were conducted using the Trimmomatic and PRINSEQ software (Bolger, Lohse & Usadel, 2014; Schmieder & Edwards, 2011). To assemble the overlapping Illumina Paired-end reads PANDAseq was used (Masella et al., 2012). In order to determine operational taxonomic units (OTUs), sequences sharing 97% identity were assembled as suggested by Kunin et al. (2010); this was done using the software QIIME (Caporaso et al., 2010). The OTUs were aligned using the GreenGenes database (http://greengenes.lbl.gov). Bacterial diversity was studied using the Shannon diversity index calculated for each DNA pool. The relative abundance of each OTU was estimated by examining the number of reads for each sequence and each sample as recommended by Jousselin et al. (2016). Taking into account that bacterial DNA contaminants can be commonly found in DNA extraction kits and other laboratory reagents or could enter samples during analysis (Salter et al., 2014), reads from taxa accounting for <1% of all the reads of a given sample were excluded from the data analysis (“unrepresented reads”). Regarding this, Jousselin et al. (2016) found that the removal of low frequency sequences (<1%) excluded the most DNA contaminants allowing for increased repeatability and reliability of results. They showed that by using this method, DNA contaminants have little impact on the analysis of aphid endosymbionts when using 16S rRNA Illumina sequencing. While, reads for which significant BLAST hits with known taxon could not be found are detailed as “unassigned reads”.

Identifying Pseudomonas species by 16S Sanger sequencing

From the 16S rRNA amplicon sequencing, a species of Pseudomonas was encountered (see ‘Results’). In order to characterize the Pseudomonas species from the 16S rRNA sequences identified, a portion of the 16S and 23S ribosomal genes (∼1,500 bp) was amplified and sequenced in 20 aphids collect from field and used to prepare sample SA-1 (Table 1); this was done using the universal bacterial primers 10F and 35R (Sandström et al., 2001; Russell & Moran, 2005). These primers were selected because they target the intergenic spacer between the 16S and 23S genes, which can be used to avoid amplifying the aphid primary endosymbiont, B. aphidicola as both genes are not contiguous in this endosymbiont (Russell & Moran, 2005).

The PCR reactions were performed in a total volume of 25 µl including; 2.5 µl of 10 × buffer, 0.2 mM dNTP’s, 2 mM MgCl₂, 0.3 µl of Taq (5 U/µl), each primer at 0.5 uM, and 3 µl of DNA (10 ng/µl). The PCR conditions consisted of initial denaturation at 94 °C for 5 min, followed by 35 cycles of 94 °C for 40 s, 57 °C for 40 s and 72 °C for 45 min and a final extension at 72 °C for 7 min. The resulting amplicons were sequenced in an ABI PRISM® 310 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). The alignment of the sequences with known 16S rRNA sequences from species of the genus Pseudomonas was conducted in Geneious v.8.1 (Drummond et al., 2011). All sequences of the genus Pseudomonas were obtained from GenBank, including sequences from the seven Pseudomonas clusters reported by Anzai et al. (2000): “Pseudomonas syringae group”, “Pseudomonas chlororaphis group”, “Pseudomonas fluorescens group”, “Pseudomonas putida group”, “Pseudomonas stutzeri group”, “Pseudomonas aeruginosa group”, and “Pseudomonas pertucinogena group” (Data S1). A phylogenetic tree of the Pseudomonas sequences was constructed using the HKY genetic distance model and the neighbor-joining method implemented in Geneious v.8.11 (Drummond et al., 2011). Branch significance was calculated using bootstrap values of 1,000 replications.

16S rRNA amplicon sequencing

A total of 1,327,786 reads were obtained after filtering the four DNA pools of the aphid S. avenae (SA-1, SA-2, SA-3 and SA-4) (Table 1). The mean Shannon diversity index was 1.51 (SD = 0.61) and ranged from 1.13 to 2.39 for the aphid S. avenae (Table 1). Of the total reads for S. avenae, 98% were classified as Gammaproteobacteria and included mostly bacteria from the Enterobacteriaceae (94.7% of the total reads) (Buchnera aphidicola, Regiella insecticola and Hamiltonella defensa) and to a lesser extent from the Pseudomonadaceae (Pseudomonas) families (3.3% of the total reads). The Gammaproteobacteria, the aphid primary endosymbiont, Buchnera aphidicola, was the most common endosymbiont in the four samples (84.4% of the total reads). The second most common taxon was the aphid secondary endosymbiont R. insecticola (9.3% of the total reads), which was also found in all studied samples (Fig. 1). Whilst another well-known aphid endosymbiont, H. defensa, represented an average of 0.9% of the total reads in two of the four samples studied (SA-1 and SA-3) (Fig. 1 and Data S2). Pseudomonas sp. sequences were well represented in two aphid samples (SA-1 and SA-4), making up an average of 3.3% of the total reads (Fig. 1 and Data S2). Unrepresented reads (i.e., reads of taxa accounting for <1% of all the reads; see methods) were found in an average of 1.5% of the total reads. Also, the four DNA pools of S. avenae had a low proportion of unassigned reads (i.e., reads for which no significant BLAST hits with known taxon were found); unassigned reads ranged from 0.5% to 0.6% (average of 0.5% of the total of reads).

For the aphid R. padi, a total of 2,095,602 reads were obtained from the six DNA pools analyzed (RP-1, RP-2, RP-3, RP-4, RP-5 and RP-6) (Table 1). A lower bacterial diversity than S. avenae was observed with an estimated mean Shannon diversity index of 0.07 (SD = 0.04) (Table 1). B. aphidicola was found in a percentage >98.5% in all DNA pools, however no additional bacteria were found in R. padi (Fig. 1). Finally, a low proportion of unassigned reads was detected among the six DNA pools (<0.01%), as well the proportion of unrepresented reads was found in an average of 0.24% of the total reads, being the highest proportion of unrepresented reads detected in the DNA pool RP-1 (1.3% of the total of reads) (Fig. 1). Sequencing data generated on Illumina were submitted to GenBank (accession numbers: MG958610, MG958611, MG958612 MG958613 and MG958614) (Data S2).

16S rRNA sequencing and phylogenetic analysis of Pseudomonas species

Of the sequences generated for the 20 aphid samples of S. avenae, only one DNA sample corresponded to a Pseudomonas species (GenBank accession number MF536106). Sequences of the other DNA samples were observed as belonging to some of the other aphid secondary endosymbionts (R. insecticola and H. defensa), as it was identified by the 16S rRNA sequencing. The phylogenetic tree constructed show the seven clusters previously described for the genus Pseudomonas (Fig. 2). The Pseudomonas sequence generated from aphid DNA was located into the “P. fluorescens group”, being closely related to Pseudomonas palleroniana with an identity percentage >95% (Fig. 2).

Secondary endosymbionts in the introduced aphid populations

A low bacterial diversity in the introduced populations of the cereal aphids S. avenae and R. padi was revealed by 16S rRNA amplicon sequencing in Chile. Gammaproteobacteria was the most common class identified and as expected the aphid primary endosymbiont, B. aphidicola, was the most common bacterial species detected in S. avenae and R. padi. In all DNA pools of both aphid species, Buchnera made up a large percentage of all of the reads (ranged between 84.4% and 99% respectively). In contrast to our systems, a greater diversity of secondary endosymbionts can be found in other aphid species (Zytynska & Weisser, 2016). For instance, the well-studied pea aphid, A. pisum, hosts at least eight secondary endosymbionts (Serratia symbiotica, R. insecticola, H. defensa, Rickettsiella, PAXS, Spiroplasma, Rickettsia and Wolbachia) that are highly abundant according to two 16 rRNA amplicon sequencing studies (Russell et al., 2013; Gauthier et al., 2015). A similar study in several aphid species of the genus Cinara has detected the presence of four dominant secondary endosymbionts (S. symbiotica, H. defensa, R. insecticola and Wolbachia) (Jousselin et al., 2016). Our results show an absence of secondary endosymbionts in the introduced populations of the aphid R. padi. Whilst the diversity of secondary endosymbionts observed in the aphid S. avenae was low; only two well-known aphid endosymbionts were encountered, R. insecticola and H. defensa. A similar bacterial diversity has been found in other geographic populations of these aphid species. For instance, aphid samples screened for a set of endosymbionts (H. defensa, R. insecticola, S. symbiotica, Rickettsia, Rickettsiella and Spiroplasma) have found only three secondary endosymbionts (H. defensa, R. insecticola and S. symbiotica) in the native range of S. avenae (U.K. and Germany) (Łukasik et al., 2013; Henry et al., 2015; Alkhedir et al., 2015). In particular, a positive association between H. defensa and S. avenae was found, being this the most common endosymbiont followed by R. insecticola, whilst S. symbiotica was reported in a lower frequency (≤6%) in the aphid populations (Łukasik et al., 2013; Henry et al., 2015). Differently, higher infection rates of R. insecticola and S. symbiotica were found in Chinese populations of S. avenae (Luo et al., 2016), as well a high prevalence of R. insecticola (75% of infected aphids) was found in introduced populations of S. avenae in Morocco (Fakhour et al., 2018). In this study, we found that R. insecticola was the dominant secondary endosymbiont in S. avenae, while H. defensa was observed at lower prevalence among DNA samples studied. However, the read abundance should be interpreted carefully when it is used as an estimate of the infection frequency of endosymbionts, because PCR amplification bias can be introduced by primer specificity (Klindworth et al., 2012). Despite this, our results from the deep sequencing of 16S rRNA gene are consistent with previous PCR-based studies on Chilean populations of S. avenae, in which ∼50% of the aphids harbored R. insecticola and a lower proportion of aphids harbored H. defensa (between 4% and 15%) (Sepúlveda et al., 2016; Zepeda-Paulo, Villegas & Lavandero, 2017), suggesting that the aphid-endosymbiont associations can vary across geographic range of aphid species.

Secondary endosymbionts make up an important component of the bacterial community of aphids and several studies have indicated that they have important effects on the host phenotype. Specifically, aphid secondary endosymbionts can protect the host from natural enemies, can provide tolerance to heat shock and can facilitate the colonization of new host plants (Oliver et al., 2010). Although recent studies have not found evidence that the endosymbionts R. insecticola nor H. defensa can confer defense against parasitoid wasps in S. avenae (Łukasik et al., 2013; Zepeda-Paulo, Villegas & Lavandero, 2017), at least one strain of R. insecticola has shown to provide protection to S. avenae against the pathogenic fungus Pandora neoaphidis (Łukasik et al., 2015). This symbiont-mediated advantage could explain the higher prevalence of R. insecticola in the populations of S. avenae here studied; however, this is not consistent with the lower prevalence of this endosymbiont reported in native regions of S. avenae. An explanation for this observation could be the founder effect and drift experienced by aphid populations introduced in a new region (Desneux et al., 2018). During the invasive process only a subset of symbiont-harboring aphid clones may have been introduced from the native regions, resulting in particular aphid—endosymbiont associations in the novel established populations. Indeed, variation in the associations between aphid clones and endosymbionts can be found in field populations, suggesting that they could be relevant for understanding of aphid—symbiont populations dynamics (Zepeda-Paulo, Villegas & Lavandero, 2017). In addition, we cannot rule out the effect of sampling method (e.g., number and distribution of sampling in a season) on the infection rates observed in aphid populations, since the frequency of endosymbionts can increase and/or fluctuate during the course of a season (Henry et al., 2015). In this regard, our aphid sampling would be considered representative of the endosymbiont diversity, as it was performed during the period of highest abundance of aphids (Raymond, Ortiz-Martínez & Lavandero, 2015; Ortiz-Martínez & Lavandero, 2018) and endosymbionts of the populations in the field (F Zepeda-Paulo & B Lavandero, 2018, unpublished data).

Unlike S. avenae, there is little knowledge on the diversity of bacterial endosymbionts in the aphid R. padi. Despite this, the existing data are consistent with our results in show an absence of secondary endosymbionts in aphid samples from the native range of R. padi (Europe) analyzed using species-specific primers developed for three aphid endosymbionts (H. defensa, R. insecticola and S. symbiotica) (Henry et al., 2015; Desneux et al., 2018) nor their introduced range (Morocco) using 16S rRNA gene sequencing (Fakhour et al., 2018). The bacterial diversity could be non-randomly distributed throughout host species. In this sense, it has been raised that the prevalence of secondary endosymbionts in a particular insect host may depend on the balance between the costs and benefits of harboring symbionts (Oliver, Smith & Russell, 2014). Indeed, the lack of an important protective phenotype providing direct benefits, fitness costs on symbiont-harboring host and the transmission rates of endosymbionts are some of the factors that could explain the low occurrence of endosymbionts in a particular host species (Oliver, Smith & Russell, 2014; Dykstra et al., 2014). Another factor that may influence bacterial diversity of aphids are the symbiont-symbiont interactions, such as competition between primary and secondary endosymbionts. Regarding this, several studies have shown that the density of the aphid primary endosymbiont, B. aphidicola, could be affected by the coexistence with secondary endosymbionts in the same host (Koga, Tsuchida & Fukatsu, 2003; Sakurai et al., 2005; Leclair et al., 2017). A negative effect on Buchnera abundance may be detrimental to the fitness of aphids and could significantly affect some aphid species. Aphids species can vary in their ability to increase the amino acid concentration in the phloem, in response to chlorotic damage induced by them (Sandström, Telang & Moran, 2000). This increase may reduce the nutritional dependence of aphids on Buchnera for the synthesis of essential amino acids, which could affect the aphid-symbiont associations. For instance, R. padi could show a high dependence on Buchnera for the synthesis of essential amino acids, since this does not affect the phloem composition of the host plant, compared to a higher amino acid concentration induced by other aphid species (Sandström & Moran, 1999; Sandström, Telang & Moran, 2000). A greater dependence in Buchnera could limit the infection of secondary endosymbionts, if they affect the abundance of the primary endosymbiont of hosts and thus explain the absence of secondary endosymbionts in some aphid species. However, the association between Buchnera-dependent aphids and the prevalence of secondary endosymbionts still have to be studied for a better understanding of the role of symbiont–symbiont interactions on the bacterial diversity of aphid species.

Presence of Pseudomonas sp. in cereal aphids

In addition to the most common aphid endosymbionts, the results from 16S rRNA sequencing showed the occurrence of Pseudomonas sp. in two DNA pools analyzed of the aphid S. avenae. However, of the sequences generated for 20 aphid samples of S. avenae, only one DNA sample corresponded to a Pseudomonas species. The phylogenetic analysis incorporating known Pseudomonas sequences showed clustering with the “P. fluorescens group”; the Pseudomonas sp. sequence generated here was closely related to the bacteria P. palleroniana and P. tolassi. These bacterial species are known phytopathogenic Pseudomonads, which have been found in rice (Oryza sativa) and garlic (Allium sativum), respectively (Gardan et al., 2002; Höfte & De Vos, 2007). Others studies based on 16 rRNA amplicon sequencing have identified phytopathogenic Pseudomonas sp. in the pea aphid (Pseudomonas syringae) and R. padi (P. viridiflava and P. veronii) (Gauthier et al., 2015). Moreover, the pea aphid has previously proven capable of acting as both a vector and a non-plant host for P. syringae (Stavrinides, McCloskey & Ochman, 2009). Some strains of P. syringae could be pathogenic to aphids causing death by bacterial sepsis (Stavrinides, McCloskey & Ochman, 2009; Hendry, Clark & Baltrus, 2016). The finding of Pseudomonas sp. in different aphid species suggests that these types of phytopathogen-vector associations may be more common than previously thought among aphid species. Secondary endosymbionts can also influence the interactions between phytopathogens and insects. Hendry, Hunter & Baltrus (2014) reported that secondary endosymbionts can influence interactions between whiteflies and the phytopathogen (P. syringae); whiteflies harboring Rickettsia have decreased their mortality from P. syringae (Hendry, Hunter & Baltrus, 2014). This latter finding might suggest that similar interactions among endosymbiotic and phytopathogenic bacteria may also occur in other host insects (Gonzalez et al., 2016). However, there are currently no studies on the extent of phytopathogen-vector/host associations or the effect of secondary endosymbionts on the interactions between aphids and phytopathogenic bacteria.