Bacterial communities of the cotton aphid Aphis gossypii associated with Bt cotton in northern China

Aphids are infected with a wide variety of endosymbionts that can confer ecologically relevant traits. However, the bacterial communities of most aphid species are still poorly characterized. This study investigated the bacterial diversity of the cotton aphid Aphis gossypii associated with Bt cotton in northern China by targeting the V4 region of the 16S rDNA using the Illumina MiSeq platform. Our sequencing data revealed that bacterial communities of A. gossypii were generally dominated by the primary symbiont Buchnera, together with the facultative symbionts Arsenophonus and Hamiltonella. To our knowledge, this is the first report documenting the facultative symbiont Hamiltonella in A. gossypii. Moreover, the bacterial community structure was similar within aphids from the same province, but distinct among those from different provinces. The taxonomic diversity of the bacterial community is greater in Hebei Province compared with in samples from Henan and Shandong Provinces. The selection pressure exerted by the different geographical locations could explain the differences found among the various provinces. These findings broaden our understanding of the interactions among aphids, endosymbionts and their environments, and provide clues to develop potential biocontrol techniques against this cotton aphid.


Comparisons of bacterial communities from different provinces. The samples from Hebei Province
were richer, having a higher number of operational taxonomic units (OTUs) than the samples from Henan and Shandong Provinces (Table 1). At the family level, the samples from Hebei generally also had more OTUs than those from Henan and Shandong, when the relative abundances of the top 35 OTUs were compared (Fig. 1). Additionally, the samples from Hebei had generally higher Ace and Chao1 richness estimates compared with the samples from Henan and Shandong (Table 1). Shannon and Simpson diversity indices also suggested that the taxonomic diversity of the bacterial community is higher in Hebei than in Henan and Shandong (Table 1). Principal coordinate analysis (PCoA) showed a distinct clustering among the individual samples, and the samples from same province tended to cluster together (Fig. 2).

Discussion
Our sequencing data revealed that the bacterial communities in A. gossypii were generally dominated by the universally present primary symbiont Buchnera. Two facultative symbionts, Arsenophonus and Hamiltonella, were also found in all of the aphid samples. To our knowledge, this is the first report finding the facultative symbiont Hamiltonella in A. gossypii. Moreover, the bacterial community structure was similar within the same province, but distinct among different provinces (Fig. 2). Our results suggest that the bacterial diversity of A. gossypii is related to the geographical location.
To utilize phloem sap as their sole dietary component, most aphids are critically dependent on symbiosis with the bacteria B. aphidicola 6 . Buchnera was also the predominant genus found in the bacterial communities of our aphid samples. Previous studies have shown that facultative symbionts from five genera can infect A. gossypii: Arsenophonus, Regiella, Rickettsia, Serratia and Wolbachia [18][19][20][21][22][23][24] . Among these, Arsenophonus was found in all of our samples, Wolbachia was only detected in some of our samples, and they both had low relative abundances. The other three symbionts were not detected in our sequencing data. Brady et al. summarized the infection rates of Regiella, Rickettsia and Serratia in A. gossypii, and each of them was only 1% 25 . This may explain why these three facultative symbionts were not found in our results.
Arsenophonus is widespread in numerous insect species 26 . It can manipulate the reproduction of various parasitoid wasps by distorting the progeny sex ratio towards the production of females through male killing 26,27 . Arsenophonus may also behave as an obligate mutualist in hematophagous insects 28 , or as a facultative mutualist, protecting against parasitoid attacks in psyllids 29 . Jousselin et al. identified aphids as harbouring an important diversity of Arsenophonus strains, and the incidence was especially high in the Aphis genus 24 . Jousselin et al. also speculate that plant mediation and parasitism might be involved in the dispersal of Arsenophonus 24 . Moreover, Arsenophonus was reported to be involved in host plant specialization in the polyphagous aphid, Aphis craccivora 30 . In recent studies, Arsenophonus did not defend its aphid host Aphis glycines against major parasitoids and fungal natural enemies 31 , but provided a general benefit to A. glycines 32 .
Wolbachia is typically associated with manipulating the reproduction of several arthropod hosts, and it can rapidly reach a high frequency in host species as a consequence 14,33 . Additionally, Wolbachia has been implicated in providing a defence against viruses in other insects 34 , and it also can provide vitamin B to the host insect 35 . In our study, Wolbachia was only found in some of the aphid samples and their relative abundance was extremely low. Similarly, Liu et al. found a low titre of Wolbachia in A. glycines 36 . The detection of Wolbachia in aphids has some difficulties, which likely has resulted in it being under reported. One major difficulty is the current PCR protocols for the detection of Wolbachia were far from optimal 18 . In addition, the development of efficient Wolbachia detection was hindered by the presence of low titre infections and multiple infections 37,38 . Wolbachia density may be affected by co-infection with other Wolbachia strains or other vertically transmitted symbionts, as well as by host genotype 39,40 . Another difficulty in detection of Wolbachia is horizontal transfer of Wolbachia genes to host genomes 41,42 , which further complicates Wolbachia detection.
Surprisingly, Hamiltonella, which has not been reported in A. gossypii, was found in our data. Hamiltonella is known to protect aphids against parasitism. Oliver et al. found that Hamiltonella defensa reduced the rate of successful parasitism by Aphidius ervi by killing developing wasp larvae 43 . Multiple strains of H. defensa were examined in A. pisum and all of them conferred a partial protection against attack by A. ervi, indicating that H. defensa generally provides protection against this wasp 10 . In addition, parasitized A. pisum containing H. defensa produced significantly more offspring than parasitized uninfected aphids, indicating that A. pisum benefited from the H. defensa infection when under attack by parasitoids 44 . H. defensa being present in all of our samples may also suggest the importance of its role in protecting A. gossypii against parasitism in the field, and both vertical and horizontal transmission may act as drivers of Hamiltonella dispersal.
Many other bacterial taxa were also detected in some of the aphid samples, and their relative abundances were extremely low. For example, Stenotrophomonas, Brevundimonas and Burkholderia, which are commonly detected in environmental samples, were also found in our aphid samples 22,23 . Burkholderia is present in the environment, associated with insects and, in some instances, clearly acts as a mutualist 45 . Many of these bacteria could be contaminants. Additionally, our aphid samples from 10 different field locations could encounter different environmental factors. Although heritability cannot be ruled out, it is more likely that these bacteria engage in opportunistic associations with aphids (perhaps as gut associates or pathogens) or that they represent contaminants from soil, plants or human handling.
Geographical location and environmental factors may account for the different bacterial community structures found in aphids from different provinces. Natural populations of aphids may experience selection pressures from various management practices, natural enemies (pathogens, predators and parasitoids) and environmental conditions that could alter the composition, as well as the frequency, of associated bacteria 46 . In the PCoA, samples from the same province tended to cluster together (Fig. 2), suggesting that the samples might have experienced similar selection pressures.
In our results, the samples from Hebei Province had a higher number of OTUs, higher Ace and Chao1 richness estimates, and higher Shannon and Simpson diversity indices, compared with the samples from Henan and Shandong Provinces (Table 1). This may be because Hebei Province has a relatively higher latitude and may encounter different climactic variables, including temperature and precipitation. A geographical-based variation in infection frequencies has been reported for some facultative symbionts. In the pea aphid, A. pisum, symbiont prevalence was found to correlate with climactic variables 47 . Bacteria from the whitefly Bemisia tabaci differ in frequency based on the host plants and geographical locations 48 . Moreover, geographical variation in Arsenophonus symbiont prevalence was reported in the psyllid Glycaspis brimblecombei 29 . These studies suggest that environmental factors may act as important drivers of natural symbiont dynamics.
Here, we used the deep Illumina MiSeq sequencing of 16S rDNA genes, to analyse the bacterial communities of the cotton aphid A. gossypii associated with Bt cotton in northern China. Our sequencing data revealed that the bacterial communities of A. gossypii were dominated by the primary symbiont Buchnera, together with facultative symbionts that varied in incidence among the aphid samples. To our knowledge, this is the first report documenting the facultative symbiont Hamiltonella in A. gossypii. Selective pressures exerted by the geographical location could explain why the bacterial community structure was similar within the same province, but distinct among different provinces. These findings increase our understanding of the intricate symbiotic relationships between bacteria and A. gossypii. Further studies will be focused on identifying the functions of the representative bacterial species and determining whether these species could play important roles in the future as biocontrol agents.

Methods
Insect sampling and DNA extraction. Apterous adults of Aphis gossypii were collected from 10 field populations in northern China during August 2014 (Table 2 and Fig. 3). All of the sampling sites were planted with Cry1Ac cotton. The aphids were immediately immersed in 90% ethanol and frozen at −80 °C upon return to the laboratory.
Prior to DNA extractions, aphid samples, each comprising 20 adult aphids, were washed for 5 min in 70% ethanol and rinsed three times with sterile water to remove surface contaminants. Then, samples were  products were mixed in equidensity ratios and mixture of PCR products was purified using the GeneJET Gel Extraction Kit (Thermo Scientific, USA).
Sequencing libraries were generated using a NEB Next Ultra DNA Library Prep Kit for Illumina (New England Biolabs, UK). The final quality and concentration of each library were checked using Agilent 2100 Bioanalyzer Instruments (Agilent Technologies, USA) and determined using KAPA Library Quantification Kits (Kapa Biosystems, USA). Sequencing was conducted on an Illumina MiSeq 2 × 250 platform at Novogene Bioinformatics Technology (Beijing, China) according to protocols described by Caporaso et al. 49 and Kozich et al. 50 .
Bioinformatics and statistical analysis. Paired-end reads were assigned to samples based on their unique barcodes and truncated by cutting off the barcode and primer sequence. Then, the paired-end reads were merged into single, longer sequences using FLASH (Version 1.2.7) 51 . Quality filtering on the raw tags was performed under specific filtering conditions to obtain high-quality clean tags 52 according to the QIIME (Version 1.7.0) 15 quality controlled process. Chimeric sequences were detected and removed using the UCHIME algorithm 53 .
Sequence analyses were performed using Uparse (Version 7.0.1001) 54 . Sequences with ≥ 97% similarity were assigned to the same OTU. Representative sequences from each OTU were screened for further annotation. For each representative sequence, the GreenGene Database 55 was used with the RDP classifier (Version 2.2) 56 to annotate taxonomic information. To study the phylogenetic relationships of different OTUs, and the differences in the dominant species of different samples, multiple sequence alignments were conducted using MUSCLE (Version 3.8.31) 57 .
To account for inequalities in sequence read depths among the samples, a randomly selected subset of 10,435 sequences per sample was chosen for a further bacterial community analysis. The microbial diversity was analysed using QIIME V1.7.0 and displayed with R software (Version 2.15.3) 15 . Alpha diversity analysis included observed species, Ace and Chao1 estimators, Simpson and Shannon diversity indices and Good's estimate of coverage. A PCoA 16 was performed to explore the differences in the bacterial community structures and was displayed with the WGCNA, stat and ggplot2 packages in the R software (Version 2.15.3). The sequencing data has been submitted to the NCBI database as a file under accession number SRP066541.