The Divergence in Bacterial Components Associated with Bactrocera dorsalis across Developmental Stages

Eco-evolutionary dynamics of microbiotas at the macroscale level are largely driven by ecological variables. The diet and living environment of the oriental fruit fly, Bactrocera dorsalis, diversify during development, providing a natural system to explore convergence, divergence, and repeatability in patterns of microbiota dynamics as a function of the host diet, phylogeny, and environment. Here, we characterized the microbiotas of 47 B. dorsalis individuals from three distinct populations by 16S rRNA amplicon sequencing. A significant deviation was found within the larvae, pupae, and adults of each population. Pupae were characterized by an increased bacterial taxonomic and functional diversity. Principal components analysis showed that the microbiotas of larvae, pupae, and adults clearly separated into three clusters. Acetobacteraceae, Lactobacillaceae, and Enterobacteriaceae were the predominant families in larval and adult samples, and PICRUSt analysis indicated that phosphoglycerate mutases and transketolases were significantly enriched in larvae, while phosphoglycerate mutases, transketolases, and proteases were significantly enriched in adults, which may support the digestive function of the microbiotas in larvae and adults. The abundances of Intrasporangiaceae, Dermabacteraceae (mainly Brachybacterium) and Brevibacteriaceae (mainly Brevibacterium) were significantly higher in pupae, and the antibiotic transport system ATP-binding protein and antibiotic transport system permease protein pathways were significantly enriched there as well, indicating the defensive function of microbiotas in pupae. Overall, differences in the microbiotas of the larvae, pupae, and adults are likely to contribute to differences in nutrient assimilation and living environments.


INTRODUCTION
Many microorganisms reside on the insect exoskeleton, in the gut and hemocoel, and within insect cells (Douglas, 2015), and relationships ranging from parasitism to mutualism are built between microorganisms and insects (Berasategui et al., 2016). These microorganisms are often identified as symbionts of insects (Douglas, 2015). Most of the best-described associations among mutualisms are based on nutritional or defensive services provided by the symbionts to their hosts. The host's physiology and behavior are often affected in such mutualisms, and the adaptability of the hosts is increased (Ohkuma and Brune, 2011;Ye et al., 2014). The symbionts provide nutrients, such as amino acids and vitamins, or digestive enzymes that aid in the degradation of fastidious dietary polymers or in the detoxification of noxious secondary metabolites (Douglas, 2009). The microorganisms can also protect their host against pathogens, parasites, parasitoids, or predators by producing toxins or antimicrobial compounds for defense (Teixeira et al., 2008;Florez et al., 2015;Hamby and Becher, 2016). Insects can generally employ the symbiont-produced defensive compounds in two different manners: (i) for the protection of the host or its offspring against antagonistic micro-or macroorganisms or (ii) as weed killers in insect fungiculture (Kaltenpoth, 2009;Ramadhar et al., 2014). Antimicrobial compounds produced by the defensive symbiont are of particular importance to insects living in enclosed, humid environments, where opportunistic fungal or bacterial infections can develop rapidly (Douglas, 2015). Studies have recently indicated that metabolic and adaptive abilities allow different bacteria to occupy their host during different host development stages and that the relationships can be multifaceted, varying in their impact on host biology (Lindow and Brandl, 2003;Turnbaugh et al., 2007;Knief et al., 2012). Hosts thus often exploit beneficial symbioses to augment their functional capabilities and to facilitate their adaptation to novel niches (Rio et al., 2006;Ye et al., 2014;Rafael et al., 2016).
For fruit fly, the pivotal roles of microbiota have been identified in recent years, and the factors that affected the structure of microbiota were also investigated. For example, the microbiotas during ontogeny of Ceratitis capitata have been reported to be shaped by phylogenetic, metabolic, and taxonomic diversities (Aharon et al., 2013). Some probiotics can even improve the fitness sexual performance of the males at emergence (Hamden et al., 2013). In Drosophila melanogaster, symbiotic bacteria play a role in mating preference by changing the levels of cuticular hydrocarbon sex pheromones (Sharon et al., 2010). For Bactrocera dorsalis, many studies have identified the structure and function of the gut microbioa. By culture-dependent and the 16S rRNA sequencing methods, the diversity of the cultivable gut bacterial communities associated with B. dorsalis have recently been investigated (Wang et al., 2011Gujjar et al., 2017), and the development and drug resistantance of B. dorsalis were affected by the gut symbionts (Cheng et al., 2017;Khaeso et al., 2017).
Studies have indicated that diet and environment greatly influence the structure of microbiotas (Egert et al., 2004;Antwis et al., 2014;Rebollar et al., 2016). Microbial communities from the surrounding environment can even serve as reservoirs and source pools of colonizers (Kueneman et al., 2014;Loudon et al., 2014). B. dorsalis undergoes great changes in living environment during its life, as eggs and larvae in fruit (it is of great possibility to be infected by microbes for B. dorsalis larvae in the rotten fruits), pupae (especially in the early stage of pupation) in the ground (enclosed, humid environments, where opportunistic fungal (especially the Metarhizium and Beauveria) can develop rapidly (Vänninen et al., 2000) and adults on the branches of fruit trees. In addition, larvae must feed on food with high sugar content, the adults must feed on food with high sugar and protein content, and pupae do not eat. And we even found the control efficiency of Metarhizium and Beauveria to B. dorsalis in the pupal stage is very low (data unpublished). These traits may result in differences in the microbiotas of B. dorsalis, and B. dorsalis may rely on multiple microbial species for fitness and provide a unique model to investigate and compare the population dynamics of symbionts that display varying levels of integration with host biology. Available data on the microbiota of B. dorsalis are unfortunately limited, which restricts understanding the microbiota's influence on host traits, such as diet and living environment.
The larvae and adults of B. dorsalis must feed on large amounts of high sugar content food, and the larvae and pupae are exposed to environments with many pathogenic microorganisms. We thus proposed the hypothesis that the symbionts of B. dorsalis will change with the development stages: in the larval and adult stages, symbionts promote the host's absorption of sugar, and some symbionts in the larval and pupal stages may also generate antibiotics to enhance resistance to the pathogens. We explored three questions in the current study. We first examined whether differences in the living environments between adults, larvae, and pupae of B. dorsalis are reflected in differences in their bacterial communities, for example more defensive bacteria in larvae and pupae. Second, as larvae and adults must feed on high sugar content food, we tested the hypothesis that functional gene abundances in larval and adult microbiotas would reflect their capacities for sugar and protein metabolism. Finally, we investigated whether functional gene abundances in larval and pupal microbiotas reflect the resistance to pathogens.

Rearing and Collection of B. dorsalis
The lab population of B. dorsalis was collected from a carambola (Averrhoa carambola) orchard (N 23 • 06 53.09 , E 113 • 24 51.29 ) in Guangzhou, Guangdong Province in April 2008 and was reared as previously described (Cheng et al., 2017). Briefly, the flies were reared under the following conditions: 25 ± 1 • C; 16:8 h light:dark cycle; 70-80% relative humidity (RH). The flies were reared with artificial diets which were treated with high pressure sterilization. Larvae, pupae, and adults of two wild populations were also collected from the cities Huizhou (HZ population) (N 23 • 25 56.00 , E 114 • , 28 16.61 ) and Nansha (NS population) (N 22 • 42 25.81 , E 113 • 33 6.30 ) of Guangdong Province in June 2017. For wild populations, carambolas with larvae were collected and take into the lab. Then pupation and eclosion processes went on in the sterile sands. Seventeen samples were collected for the lab population (six larvae, six pupae, and five adults); 15 samples were collected for the HZ and NS population (five larvae, five pupae, and five adults). Each sample consists of one individual.

Bacterial Community Characterization
All samples (Whole individuals) were selected and then washed with sterile water. The washed samples were transferred to centrifuge tubes containing DNA extraction buffer (with lysozyme) and grinded with a homogenizer. Then total DNA of the samples was extracted using a DNA extraction kit (Tiangen, Beijing, China) following the manufacturer's instructions. After the DNA of the samples was prepared, qPCR was used to estimate the bacteria absolute content of the samples with the universal bacterial 16S rRNA primers. A standard curve for qPCR was generated by amplifying the 16S rDNA of the Arthrobacter sp. isolated from the pupa of B. dorsalis. Approximately 465 bp of the V3-V4 region of the bacterial 16S rDNA gene was amplified by PCR according to a standard protocol. The following primers were used: F, 5 -CCTACGGGNGGCWGCAG-3 and R, 5 -GGACTACHVGGGTATCTAAT-3 . The primers contained the A and B adapters for 454 Life Sciences pyrosequencing and a unique 12-bp error-correcting Golay barcode, which allowed multiplexing of samples in a single run. Each sample was analyzed in a total reaction volume of 25 µL that contained 2.5 µL of Takara 10× Ex Taq buffer, 1.5 µL of Mg 2+ (25 mM), 2 µL of dNTP (2.5 mM), 0.25 µL of Takara Ex Taq (2.5 U/µL), 0.5 µL of each primer (10 µM), 16.75 µL of ddH 2 O, and 1 µL of template. The PCR amplifications were performed with a 2-min incubation at 95 • C followed by 30 cycles of 94 • C for 30 s, 57 • C for 30 s and 72 • C for 30 s, and a final 5-min extension at 72 • C. The PCR products were purified using QIAGEN MinElute PCR Purification Kit (QIAGEN, Hilden, Germany) to remove unincorporated primers and nucleotides. A microspectrophotometer ND-1000 (NanoDrop Technologies, Wilmington, DE, United States) was used to measure the concentration of the purified DNA. The purified DNA was sequenced using an Illumina sequencing kit and an Illumina MiSeq sequencer (Illumina, San Diego, CA, United States).
Paired Illumina reads were merged in QIIME (Caporaso et al., 2010). After the high-quality reads were obtained, the data were filtered to remove low-complexity sequences (such as poly-A sequences) and sequences with ambiguous nucleotides, and the operational taxonomic units (OTUs) were picked using USEARCH (Edgar, 2010). The number of OTUs was calculated with mothur software (Schloss et al., 2009) at 97% similarity. An RDP classifier (Huse et al., 2010) was used with naïve Bayes settings for species annotation; the confidence threshold was set to 0.5. Using the species annotations and the read numbers of the OTUs, we generated OTU abundance profiles for all samples. OTUs with an abundance <0.005% of the total data set were removed as an additional level of quality filtering (Bokulich et al., 2013;Navas-Molina et al., 2013).

Diversity Analyses
Alpha diversity and Shannon rarefaction curves were calculated for all samples in mothur to investigate the species richness of the samples (v.1.34.0) (Schloss et al., 2009). Bray-Curtis and unweighted UniFrac distance matrices were used to calculate the beta diversity and were visualized with principal components analysis (PCA). To determine whether bacterial communities differed among host species, we conducted a shared and unique OTU analysis on the basis of an OTU table generated by QIIME. We used the unweighted pair group method with arithmetic mean (UPGMA), a hierarchical clustering method based on the arithmetic mean, to determine clustering patterns across host species. The UPGMA was used on Bray-Curtis distances of mean OTU relative abundances at the family level. The UPGMA, Bray-Curtis calculations and resulting heatmap were completed using the vegan package (Oksanen et al., 2015) in the R statistical package. Putative microbiota functions were predicted by annotating pathways of OTUs against the KEGG database using PICRUSt (Langille et al., 2013).

Statistical Analysis
The differences between treatments were compared by a oneway analysis of variance (ANOVA), followed by Tukey's test for multiple comparisons. The differences were considered significant at the P < 0.05 level. The data were analyzed using SPSS software. Analysis of similarity (ANOSIM) for the bacterial community of B. dorsalis across developmental stages were done with PRIMER 7.0.

Symbionts Content Quantification and Sequencing Data of 16S rRNA
Absolute content of the symbionts in flies were identified with qPCR, the results showed the symbionts content in each individual was about 10 6 CFU and there is no difference for the absolute content of the symbionts in different development stage (lab: F 2,14 = 0.126, P = 0.883; Huizhou: F 2,12 = 0.717, P = 0.508; and Nansha: F 2,12 = 1.768, P = 0.212) (Supplementary Figure S1). After the sequencing data were subjected to demultiplexing, quality filtering and chimera removal, 50082-57433 reads were retained for the 17 lab population samples, 50176-57172 reads were obtained for the 15 YC samples, and 50037-56599 reads were obtained for the 15 NS samples (Supplementary Table S1). Shannon rarefaction curves for all samples showed a plateau stage, indicating adequate sampling of 16S rRNA sequences for all the samples (Supplementary Figure S2).

Differences in Larval, Pupal, and Adult Bacterial Communities
Significantly more OTUs were generated in the pupal samples of the three populations (lab: F 2,14 = 30.387, P < 0.001; Huizhou: F 2,12 = 5.746, P = 0.018; and Nansha: F 2,12 = 5.116, P = 0.025) (Supplementary Figure S3). The ACE value and the Shannon and Simpson indices indicated that pupae exhibited the greatest species richness and that the richness of adults and larvae did not significantly differ (Supplementary Table S2). A major trend is clearly seen during development stages: the bacterial diversities of all populations were closest, on average, during the feeding stage (larvae and adults), while pupae in complex environments showed greater bacterial diversity (P < 0.01).
Larvae and adults were significantly enriched in Proteobacteria at the phylum level more than pupae (for larvae: a cumulative relative abundance above 60, 44.3, and 91.9% in larvae for the lab, Huizhou and Nansha populations, respectively, P < 0.01; and for adults: a cumulative relative abundance above 57.82, 54.32, and 92.12% for the lab, Huizhou and Nansha populations, respectively, P < 0.01). Actinobacteria had a cumulative relative abundance in pupae above 20.54, 26.01, and 19.23% for the lab, Huizhou and Nansha populations, respectively (Figure 1).
The abundance of the dominant OTUs in each stage was compared with that in the other two stages (Figure 2 and Supplementary Data Sheet S1).
In conclusion, we found Acetobacteraceae (Acetobacter sp.) was the most abundant OTU in larvae of the three populations; Intrasporangiaceae, Dermabacteraceae (mainly Brachybacterium), and Brevibacteriaceae (mainly Brevibacterium) were the most abundant OTUs in pupa of the three populations; Enterobacteriaceae was the most abundant OTU in adult of the three populations (Figure 2 and Supplementary Data Sheet S1).
Bacterial communities of larvae, pupae, and adults showed a clear pattern of specialization based on unweighted UniFrac distances with OTUs annotated at the genus level (PCA, Figure 3), indicating the specialization of larvae, pupae, and adults in hosting OTUs unique to each stage. PCA was also used to compare the similarity in the microbial community compositions of all samples of the three populations. Larvae, pupae, and adults each formed a distinct cluster among all samples, and these three clusters were separated from each other (Figure 4). The clustering pattern among samples was not influenced by the sampling population, as samples from the three populations clustered together, with the exception that results from pupal samples formed two different clusters. Moreover, ANOSIM results indicated that there were significant differences in the bacterial community of B. dorsalis across developmental stages (lab: R = 0.998, P = 0.001; Huizhou: R = 0.745, P = 0.001; and Nansha: R = 0.951, P = 0.001).

Functional Prediction of Larval, Pupal, and Adult Microbiotas
We addressed whether increased OTU diversity confers the host with a higher functional diversity. We predicted 4364, 4590, and 4308 level 3 KEGG Orthology (KO) groups in the predicted metagenomes of the lab, Huizhou and Nansha populations, respectively (PICRUSt analysis, Supplementary Data Sheet S2). The pattern of functional diversity largely followed the trend in taxonomic diversity: pupae displayed greater bacterial taxonomic and functional diversity than larvae and adults (lab: R 2 = 0.587,  P < 0.01; Huizhou: R 2 = 0.678, P = 0.005; and Nansha: R 2 = 0.3987, P = 0.012; Pearson relationship, Figure 5). The microbial community clusters of larvae, pupae, and adults were significantly separated at both the OTU and KO levels (Figures 3, 6), indicating a strong development effect.

DISCUSSION
Although studies on fly microbiotas have been reported (Aharon et al., 2013;Aksoy et al., 2014;Augustinos et al., 2015;Michael et al., 2016;Cheng et al., 2017;Yong et al., 2017b;Zhao et al., 2017), little is known about the microbial community differences during different development stages, and the former studies mainly focused on the gut microbiotas of the flies (Aksoy et al., 2014;Augustinos et al., 2015;Michael et al., 2016;Cheng et al., 2017;Zhao et al., 2017). The number of OTUs observed in the larvae, pupae, and adults of B. dorsalis in this study is greater than those reported in other fly gut samples (Figure 1), suggesting that the function of the microbial community in B. dorsalis must be analyzed in detail. Actually microbiomes associated with B. dorsalis have been reported in earlier studies. For example, Andongma et al. (2015) have reported that the dominance of Firmicutes in adult stages and Proteobacteria in immature stages (Andongma et al., 2015); Wang et al. (2011) and Yong et al. (2017b) revealed Proteobacteria (specifically Gammaproteobacteria) to be predominant in the male adults of B. dorsalis (Wang et al., 2011;Yong et al., 2017b). In our study Proteobacteria was the dominant phylum in larvae and adults, which is consistent with previous studies (Aksoy et al., 2014;Yong et al., 2017b); however, Actinobacteria were mainly found in pupae, which may indicate that the pupal microbiota has a different function. Although Andongma et al. (2015) investigated symbiotic bacterial populations across life stages of B. dorsalis, many fewer OTUs were identified in their study; moreover, fewer Actinobacteria were not identified in pupa in their study (Andongma et al., 2015). The DNA extraction method may FIGURE 7 | Comparison of predicted KEGG ortholog counts. Means ± SEMs labeled with different letters are significantly different. OTU numbers annotated for phosphoglycerate mutase (A), antibiotic transport system ATP-binding proteins (B), antibiotic transport system permease proteins (C), transketolase (D), and protease (E).
be the key reason for differences between their study and the current study. DNA is difficult to extract from Actinobacteria without lysozyme digestion because of its special cell wall (Zhang et al., 2013). Andongma et al. (2015) did not digest the sample with lysozyme, possibly causing the absence of Actinobacteria contributions. Moreover, they used pupae without puparium, which may explain the lack of Actinobacteria since bacteria of this order are often located in specific regions of the surface of insect hosts (Kaltenpoth, 2009). The greater sequencing depth in our study might also explain the greater number of identified OTUs since we obtained many more reads than Andongma reported. Actually other study has also identified Actinobacteria (Yong et al., 2017b).
Insects show remarkable adaptations to exploit diverse nutritional resources; these adaptations are due to the wide diversity of digestive enzymes produced by the insects themselves as well as the metabolic capabilities of symbiotic microorganisms that overcome the host's nutritional limitations (Berasategui et al., 2016). The high abundance of Proteobacteria observed in larvae and adults likely supports their importance in sugar metabolism. Acetobacteraceae, Lactobacillaceae, and Enterobacteriaceae were the most dominant families within this phylum that were observed in all larval and adult samples and have been reported to function in sugar metabolism (Kersters et al., 2006;Lambert et al., 2011;Yong et al., 2017a). However, we found that Acetobacteraceae and Lactobacillaceae were the most dominant families in larvae, while Enterobacteriaceae was the most dominant family in adults. This result may suggest that the digestion process may differ in larvae and adults, resulting in changes in the microbiota composition. Unlike the larvae, the adults must also feed on a high-protein diet, and protein in their diet can even influence the mating behavior of flies (Shelly and Kennelly, 2002;Shelly et al., 2005). An important family that is associated with most fruit flies is Enterobacteriaceae; members of this group play very important roles in courtship and reproduction (Ben Ami et al., 2010). We also found that pathways involved in protein metabolism were significantly enriched in microbiotas in adults. These results indicate that microbiotas may also be involved in digesting protein in food. More evidence is needed to prove whether members of Enterobacteriaceae also play a role in protein digestion to influence the courtship and reproduction of B. dorsalis. Newell et al. (2014) also reported that some Acetobacteraceae species in the gut of Drosophila were involved in oxidative stress detoxification and encoded an efflux pump (Newell et al., 2014). Researches have suspected that Lactobacillaceae and Enterobacteriaceae contribute to digestion and protection against parasites and pathogens in insect gut (Billiet et al., 2017;Smith et al., unpublished). We thus need detailed investigation of the specific bacteria of B. dorsalis by pure culture methods.
Insect-associated microbes are just beginning to be exploited as promising sources of novel bioactive compounds (Dettner, 2011). Microbial symbionts providing chemical defense for the host against predators, parasites, parasitoids, and pathogens occur in several insect taxa, including beetles (Kellner, 2002), psyllids (Nakabachi et al., 2013), planthoppers (Fredenhagen et al., 1987), and solitary wasps (Kaltenpoth et al., 2005;Kaltenpoth, 2009;Kaltenpoth and Engl, 2014). The high abundance of Actinobacteria observed in the pupae of B. dorsalis support their importance in producing compounds with antimicrobial activity. Actinobacteria are known to be important sources suited as defensive symbionts of insects (Kaltenpoth, 2009). The number of OTUs in pupae that represent antibiotic transport system ATP-binding and antibiotic transport system permease proteins is significantly greater than that in larvae and adults, which strongly supports the defensive function of Actinobacteria in the pupae of B. dorsalis. Actinobacteria are particularly common and widespread in soil (Goodfellow and Williams, 1983) and are therefore regularly encountered by insects living in soil. The pupae of B. dorsalis consistently remained in soil until emergence. Kaltenpoth (2009) stated that three key factors probably contribute to the role of Actinobacteria as defensive exosymbionts in insects: (i) their ability to utilize a wide variety of carbon sources and to generally subsist on low levels of resources; (ii) the capacity of some taxa to form spores and thereby survive inhospitable conditions in the soil; and (iii) their ability to produce secondary metabolites with antibiotic properties (Goodfellow and Williams, 1983). The evolution of symbiotic interactions between Actinobacteria and insects might have been initiated by commensal or facultative entomopathogenic Actinobacteria that exploited the low amounts of compounds present on the cuticle or in the excretions of soil-dwelling insects. Once the bacteria became associated with an insect, antibiotic substances produced for the microbes' own protection might have also become beneficial for the host insect (Kaltenpoth, 2009). The specific inhabitation of Intrasporangiaceae, Dermabacteraceae, and Brevibacteriaceae in pupae in our study may indicate their defensive function. The defensive function of bacteria in Brevibacteriaceae has been previously revealed by pure culture methods, and researchers have identified the relevant bacteria and antibacterial compound (Ryser et al., 1994;Maisnierpatin and Richard, 1995). Brachybacterium of Dermabacteraceae were also identified to express strong antimicrobial activity (Liu et al., 2011). Pupae can therefore be used in future studies as a source from which Actinomycetes with antimicrobial activity can be isolated.

CONCLUSION
The larvae, pupae, and adults of B. dorsalis were observed to harbor distinct microbial flora. We performed a detailed investigation of the microbial flora of B. dorsalis that provides a basis for future research. Further studies to investigate the microbial composition may provide a comprehensive understanding of the differences in diet and physiological behavior among B. dorsalis. Moreover, hostspecific microbial species, for example, those from the phylum Actinobacteria, can be used to develop potential compounds with antimicrobial activity that have potential value for human application.

AVAILABILITY OF DATA AND MATERIAL
Sequence data has been deposited at NCBI under Bioproject PRJNA415228.

AUTHOR CONTRIBUTIONS
Conceived and designed the experiments: DC; analyzed the data: DC and YL; contributed reagents/materials/analysis tools: XfZ, XyZ, ZC, and ZW; wrote the paper: DC.