Similar Shift Patterns in Gut Bacterial and Fungal Communities Across the Life Stages of Bactrocera minax Larvae From Two Field Populations

Bactrocera minax (Enderlein) (Diptera: Tephritidae) is an oligophagous insect pest that damages citrus fruit, especially in China. Due to larvae living within a highly septic environment, a wide variety of microorganisms exist in the larval gut of B. minax. However, a systematic study of the intestinal microbiota of this harmful insect pest is still lacking. Here, we comprehensively investigated the larval gut microbiota of B. minax in two field populations from Zigui (developed in orange) and Danjiangkou (developed in mandarin orange). We observed a dominance of Proteobacteria and Firmicutes in these bacterial communities, and Enterobacteriaceae was the predominant family throughout the larval stage. However, most of the identified fungal sequences were annotated as being from either Ascomycota or Basidiomycota phyla. Although there was a difference in the structure of the microbial communities between the two populations, the dynamic change patterns of most of the members of the microbiota were similar across the lifespan of larvae in both populations. The relative abundances of the Acetobacteraceae, Leuconostocaceae, and Lactobacillaceae gut bacteria as well as the Pichiaceae, Sebacinaceae, and Amanitaceae fungi increased throughout development, and these microorganisms stably resided in the larval gut. Furthermore, the dynamic changes of the functions of gut bacterial communities were inferred, and there was a significant increase in carbohydrate metabolism across the lifespan of larvae in both groups. Spearman correlation analysis showed that Acetobacteraceae, Lactobacillaceae, and Leuconostocaceae displayed a positive correlation with fructose and mannose metabolism, an important pathway of carbohydrate metabolism, highlighting the potential roles of these prevalent microbial communities in host biology.


INTRODUCTION
The insect gut is an ecosystem that is accompanied by complicated interactions with microbes: numerous studies have revealed the diverse compositions and structures of the microbial communities in the insect gut (Engel and Moran, 2013;Douglas, 2015), ranging from several hundred species in termites (Hongoh et al., 2005;Brune and Dietrich, 2015), to less than a hundred species in lepidopterans (Broderick et al., 2004;Robinson et al., 2010;Pinto-Tomás et al., 2011;Paniagua et al., 2018), and to almost a complete absence in aphids (Grenier et al., 2006). Gut microbes influence hostassociated biology and ecology in diverse ways, including by degrading of recalcitrant food, providing of specific nutrients (Engel and Moran, 2013), regulating energy balance (Yamada et al., 2015), and promoting host growth and development (Shin et al., 2011). Additionally, some studies have also revealed that microbiota play roles in detoxification (Ceja-Navarro et al., 2015) and in protection from predators, parasites, and pathogens (Endt et al., 2010;Stecher and Hardt, 2011), thereby influencing mating (Sharon et al., 2010), maturation and development of the immune system (Weiss et al., 2011). Changes in the intestinal microbiota composition caused by antibiotics, a dysregulated immune system or irradiation reduce the host's ecological fitness (Ryu et al., 2008;Raymann et al., 2017;Cai et al., 2018). Indeed, the insect gut microbiota is considered an essential component of host insects, and maintenance of a healthy composition of the gut microbiota is important for host life.
Bactrocera minax is an oligophagous and univoltine species that damages citrus fruit, mainly in southern/central China, and many Southeast Asian countries (Wang and Luo, 1995;Drew et al., 2006). Adult flies live by feeding on substances such as honeydew and reach sexual maturity approximately 3 weeks after emerging. Females prefer laying eggs in young fruits and oviposit from 50 to more than 200 eggs (Xia et al., 2018). Interestingly, the eggs last for approximating 1 month before hatching, and the entire lifespan of the larvae is completed in the fruit, which in turn, causes fruits to brown and rot, leading to premature drop. And B. minax experience a diapause when overwinter in the soil as pupae (Chen et al., 2016). To control this agricultural pest, a number of studies are currently being performed to understand the biology and ecology of B. minax (Xiong et al., 2016;Wang et al., 2018). However, there is still a lack of research on the bacterial and fungal compositions of the larval gut of B. minax, as well as on the shift patterns and functions of the core microbiota.
In the present study, we investigated the dynamic changes of bacterial and fungal communities in the gut of B. minax larvae at different life stages from two field populations. Apparently, the microbial communities shifted across the lifespan of the larvae in correlation with host developmental changes. Most bacteria and fungi present in the gut displayed a similar pattern of prevalence across the different larval stages of B. minax. In addition, we inferred changes in metabolic pathways during the development of larvae, which were influenced by the shifts in bacterial communities.

B. minax Larval Sample Collection
Maggot-infested fruits were collected from Zigui and Danjiangkou in Hubei Province. Before dissection, maggot fruits were surface sterilized by immersion in 70% ethanol for 3 min and rinsed three times in sterile phosphate-buffered saline to avoid contamination from the environment. The classification standard for different instars was as follows: (1) first instar larvae, body length of less than 5 mm and a milky white color; (2) second instar larvae, body length of 5-14 mm and a creamy yellow color; and (3) third instar larvae, body length of 13-18 mm and a maize-yellow color. After collecting different development stages, we dissected the whole gut from the larvae. Before dissection, the larvae were surface sterilized by immersion in 70% ethanol for 3 min and rinsed three times in sterile phosphate-buffered saline. Four independent cohorts of larvae were dissected and used as biological replicates. Tissue samples were homogenized in an automatic sample Precellys-24 homogenizer (Shanghai Jingxin Industrial Development Co., Ltd., Shanghai, China) at 70 Hz/s for 60 s with a 10 s interval.

DNA Extraction and High-Throughput Sequencing
DNA was isolated from 30 larval guts per replicate using the E.Z.N.A TM Soil DNA kit (Omega, Norcross, GA, United States) following the manufacturer's instructions. To confirm successful DNA isolation, a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) and agarose gel electrophoresis were used to check the quantity and quality of extracted DNA, respectively. The 16S ribosomal RNA (rRNA) gene-spanning variable regions (V3-V4) and fungal internal transcribed spacer (ITS1) region were amplified using the following primer pairs: 338F/806R (5 -ACTCCTACGGGAGGCAGCA-3 and 5 -GGACTACHVGGGTWTCTAAT-3 ) and ITS5/ ITS2 (5 -GGAAGTAAAAGTCGTAACAAGG-3 and 5 -GCT GCGTTCTTCATCGATGC-3 ). Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN, United States) and the PicoGreen dsDNA Assay kit (Invitrogen, Carlsbad, CA, United States) were used to purify and quantify the PCR amplicons, respectively. After quantification of individual amplicons, equal amounts of barcoded amplicons were pooled, and paired-end 2 × 300 bp reads were obtained by sequencing on the Illumina MiSeq platform at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China). The raw sequencing data were deposited in the National Center for Biotechnology Information Sequence Read Archive with accession no. PRJNA542035. The B. minax adult gut 16s rRNA data was from Wang et al. (2014) (under accession number SRR1531158 in NCBI).

Sequence Analysis and Diversity Measures
Sequencing data were processed by the Quantitative Insights into Microbial Ecology (QIIME, v1.8.0) pipeline (Caporaso et al., 2010). After filtering low-quality sequences according to previous criteria (Gill et al., 2006), FLASH was employed to assemble the paired-end reads (Magoč and Salzberg, 2011). After chimera detection, all remaining high-quality sequences were clustered into operational taxonomic units (OTUs) with a similarity threshold of 97% sequence identity by UCLUST (Edgar, 2010). The Greengenes Database (DeSantis et al., 2006) and the Unite Database (Kõljalg et al., 2013) were used for BLAST searching with the representative sequences and further taxonomy classification.
The qualified OTU data in QIIME were used to calculate the alpha diversity indices, including the Chao1 richness estimator, the abundance-based coverage estimator (ACE) metric and the Shannon and Simpson diversity index. Regardless of the relative abundance and based only on the occurrence of OTUs, the R package "VennDiagram" was used to generate a Venn diagram to visualize shared and unique OTUs across the lifespan of larvae (Zaura et al., 2009). Beta diversity analysis was employed to assess the structural variation of microbial communities among samples based on UniFrac distance metrics (Lozupone and Knight, 2005;Lozupone et al., 2007) visualized via non-metric multidimensional scaling (NMDS) and unweighted pair-group method with arithmetic means (UPGMA) hierarchical clustering (Ramette, 2007). Phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) was performed to predict microbial functions using the high-quality sequences (Langille et al., 2013).

Statistical Analysis
Student's t-test and one-way analysis of variance by Tukey's test were used to calculate the statistical significance among different groups, and a P-value < 0.05 was representative of statistical significance. Permutational multivariate analysis of variance (PERMANOVA) was performed to assess the differences in the microbiota structure among groups (Mcardle and Anderson, 2001). The fitting curve of the microbiota, Spearman correlations and correlation heatmaps were performed with R statistical software 1 . GraphPad Prism 8.0 (GraphPad Software, La Jolla, CA, United States) and R statistical software were used to generate all the diagrams in the present study.

Composition of Gut Microbiota Across the Life Stages of B. minax Larvae
High-throughput sequencing was employed to investigate the composition of the bacterial and fungal communities that inhabit the larval gut. The sequencing data yielded a total of 930,502 reads of bacterial 16S rRNA genes and 1,058,828 reads of the ITS region of fungal rRNA genes with average lengths of 435 and 202 bp, respectively. Species accumulation curve analysis indicated that the samples were sufficient to reveal the microbial communities, as the species accumulation curves tended toward saturation (Supplementary Figures S1A, S2A). The rank-abundance curves showed that most of the reads from the bacterial and fungal communities belonged to rare species and that the distribution of species was relatively concentrated according to the breadth and smoothness of the curves (Supplementary Figures S1B, S2B).
Taxonomic analysis of gut bacteria revealed that all samples were intimately associated with diverse microbes across the lifespan of larvae ( Figure 1A). The most prevalent phyla were Proteobacteria and Firmicutes, and they were detected in all larvae of B. minax from two field populations (Supplementary Figure S3). Venn diagram analysis showed that subsets of 93 1 http://cran.r-project.org/ and 112 shared bacterial OTUs resided in the gut across the lifespan of larvae from ZG and DJK, respectively ( Figure 1B). Most of the shared OTUs belonged to Enterobacteriaceae (Supplementary Data Sheet S1). Enterobacteriaceae were the core members of gut bacteria communities in both larvae ( Figure 1C) and adults of B. minax (Supplementary Table  S1). With the development of larvae, there was a significant increase in the relative abundance of Acetobacteraceae, Leuconostocaceae, and Lactobacillaceae ( Figures 1D,G,I), while the relative abundance of bacteria of these families was different in larvae at the same instar stage from different populations, especially 3rd instar larvae (Supplementary Table  S1). By contrast, we found some bacterial families, including Brucellaceae, Chitinophagaceae, Moraxellaceae, and Bacillaceae, to be decreased or absent in both field populations during the development of larvae ( Figure 1F and Supplementary  Figures S4I-K). The relative abundance of some bacterial families, including Oxalobacteraceae, Comamonadaceae, Caulobacteraceae, Enterococcaceae, Bradyrhizobiaceae, and Sphingomonadaceae, only existed in the gut of 1st or 2nd instar from ZG, and the abundance decreased to undetectable levels in 3rd instar larvae (Supplementary Figure S4). We also found a few bacteria, such as Xanthomonadaceae, varied between the two groups, and the abundance of these bacteria sharply increased in 3rd instar larvae from ZG but not from DJK ( Figure 1E). Pseudomonadaceae consistently survived in the gut of larvae collected from ZG (Supplementary Figure S4D), while Aeromonadaceae and Streptococcaceae stably resided in the gut of larvae collected from DJK ( Figure 1H and Supplementary Figure S4E). Moreover, there was a sharp decrease in the relative abundance of other bacteria (unclassified and very low abundance bacteria) in ZG samples during larval development, but the relative abundance of other bacteria was stable and persisted in the gut of larvae collected from DJK ( Figure 1J).
We further investigated the fungal composition of larval guts collected from ZG and DJK. More than 85% of the fungal sequences were annotated as belonging to the two phyla Ascomycota and Basidiomycota, which were dominant in the fungal community (Supplementary Figure S3B). There were subsets of 173 and of 35 shared fungal OTUs in the gut across the lifespan of larvae from ZG and DJK, respectively, as revealed by Venn diagram analysis ( Figure 2B). At the family level, an increase of Pichiaceae, Sebacinaceae, and Amanitaceae occurred during the development of larvae in both groups, and these fungi stably resided in the gut of 3rd instar larvae (Figures 2A,C,D,G). In addition, Pichiaceae were present across three larval stages (Supplementary Data Sheet S1). During larval development, several fungi, such as Thermoascaceae, Plectosphaerellaceae, Aspergillaceae, Metschnikowiaceae, Lasiosphaeriaceae, Cordycipitaceae, Trichocomaceae, Nectriaceae, Chaetomiaceae, Mortierellaceae, and Russulaceae, decreased or disappeared in both groups (Figures 2E,F  In every panel, the left side represents the population collected from Zigui (ZG, n = 4), and the right side represents the population collected from Danjiangkou (DJK, n = 4). The high, low and median data values are shown in boxplots, and the lower and upper edges of every box represent the first and third quartiles, respectively. Multiple comparisons were performed in the two groups separately. Different letters indicate a significant difference between different instars in each group (p < 0.05, one-way ANOVA, Tukey post hoc test). The curve fitting was created using the "ggplot2" package in R statistical analyses. The lines show median values per region window, and the shaded area denotes the estimated 95% confidence interval.
were only detected in the gut of the 2nd or 3rd instar larval population collected from ZG (Supplementary Figures  S5A,D). There was a sharp decrease of other fungi (unclassified and very low abundance fungus) in DJK samples as larval development progressed, but the unclassified and very low abundance fungi stably resided in the gut of larvae collected from ZG ( Figure 2J).

Diversity of Gut Microbiota Across the Life Stages of B. minax Larvae
The species richness of bacteria among different life stages of larvae was analyzed using the observed species (Obs), ACE, and Chao1. Diversity was analyzed using the Shannon, Simpson and Berger-Parker index. The OTUs and two richness estimators (ACE and Chao1) showed that there were no significant differences among the different development stages of larvae in either field population (Figures 3A-C). When taking into account the diversity index (Shannon index), 1st instar larvae had the greatest species diversity in both groups, and there was a decrease in species diversity with the development of larvae ( Figure 3D). While only the Simpson and Berger-Parker estimators of larvae collected from DJK displayed the same pattern as above, there were no significant differences across the lifespan of larvae collected from ZG based on the Simpson and Berger-Parker estimator (Figures 3E,F).
Analysis of the alpha diversity of fungi among the different life stages of larvae demonstrated that during development, there was a significant decrease in species richness and diversity of fungi in the gut of larvae collected from ZG based on Obs, In every panel, the left side represents the population collected from ZG (n = 4), and the right side represents the population collected from DJK (n = 4). The high, low, and median data values are shown in boxplots, and the lower and upper edges of every box represent the first and third quartiles, respectively. Multiple comparisons were performed in the two groups separately. Different letters indicate a significant difference between different instars in each group (p < 0.05, one-way ANOVA, Tukey post hoc test). The lines show median values per region window, and the shaded area denotes the estimated 95% confidence interval. ACE, Chao1, and the Shannon, Simpson and Berger-Parker index (Figures 3G-L). First instar larvae had the greatest fungal richness and species diversity. Linear curve analysis of alpha diversity showed that there was also a decrease in species richness and diversity of fungi as larvae collected from DJK developed from 1st instar to 3rd instar (Figures 3G-L), but this decrease was not significant.

Structure of Gut Microbiota Across the Life Stages of B. minax Larvae
Based on UniFrac metrics, the bacterial, and fungal community dissimilarity was investigated. UPGMA and heatmap analyses showed that there was an apparent separation between early instar (1st instar larvae) and late instar (3rd instar larvae) larvae in both field populations, and larvae from the same instar were tightly clustered (Figures 4A,B). Although 2nd instar larvae clustered with the 1st and 3rd instar larval groups, there was a significant difference in the UniFrac distance between 1st instar and 2nd instar larvae (Figures 4D,F) but no difference between 2nd instar and 3rd instar larvae (Figures 4E,G), indicating that the structure of 2nd instar larval gut bacteria is similar to that of 3rd instar larvae. NMDS analyses also revealed that the composition and structure of the gut bacterial community varied among the different life stages of larvae in both field populations. In addition, NMDS was sufficient to separate larvae of different stages and the same instar larval stage from different field populations (PERMANOVA test with 999 permutations, P < 0.001) (Figures 4C-E,H). We next analyzed fungal community dissimilarity across the lifespan of larvae. UPGMA and heatmap analyses showed that there was a relatively tight clustering of the fungal communities of larvae at the same instar stage ( Figure 5A) and that the UniFrac distances among different larval stages were significantly different in the ZG population (Figures 5D,E). There was also an apparent separation between early instar and late instar larvae in the DJK population ( Figure 5B), with a significant difference in UniFrac distance between 1st instar and 3rd instar larvae ( Figure 5F). While the 2nd instar larvae clustered with 1st and 3rd instar larval groups, there was no significant difference in the UniFrac distance between 2nd instar larvae and other development stages (Figures 5F,G). NMDS analyses also revealed that the composition and structure of the gut mycobiota varied among the different life stages of larvae in both field populations and was sufficient to differentiate the various larval stages and same instar larval stage from different populations (PERMANOVA test with 999 permutations, P < 0.001) but not from 2nd instar larvae in DJK (Figures 5C-E,H).

PICRUSt Metagenomic Predictions of Microbial Communities in Larval Gut
Approximately 41 KEGG pathways were obtained using PICRUSt to infer the functions of the bacterial communities (Supplementary Data Sheet S2). The nearest sequenced taxon index (NSTI) of larval samples is shown in Supplementary Table  S2. Here, we mainly focused on the metabolic pathways that had a >1% relative abundance. The dynamic change patterns of most metabolic pathways across the development stages of larvae were different in the two field populations (Supplementary Figure S6) except for carbohydrate metabolism, which was increased across the lifespan of larvae in both groups ( Figure 6A). Within carbohydrate metabolism, we found that several pathways were also increased with the development of larvae, such as fructose and mannose metabolism, pentose phosphate pathway, and glycolysis or gluconeogenesis (Figures 6B-D). The relative abundances of these metabolic pathways in larvae from DJK were higher than in the corresponding instar larvae from ZG (Figures 6B-D). To clarify the association between functionality and the changes in the composition of the bacterial microbiota, Spearman correlation was employed. The results showed that the Acetobacteraceae, Lactobacillaceae, Leuconostocaceae, and Xanthomonadaceae had positive correlations with almost all the predicted pathways of carbohydrate metabolism in ZG samples and that Lactobacillaceae and Leuconostocaceae displayed positive correlations with butanoate metabolism, fructose and mannose metabolism and inositol phosphate metabolism in DJK samples ( Figure 6E).

DISCUSSION
Here, we surveyed the microbiota of B. minax larvae from two different field populations, and all the samples were intimately associated with a wide variety of microbes. Similar to other In each panel, the first entry in the X axis serves as the reference for the comparisons, and the median is plotted as the horizontal line. * p < 0.05; * * p < 0.01; and * * * p < 0.001; ns, not significant as calculated by Student's t-test.
insects, we found that Proteobacteria and Firmicutes dominated the bacterial community across the lifespan of larvae (Yun et al., 2014;Chen et al., 2018). For the first time, we also investigated the mycobiota of the B. minax larval gut, and Ascomycota and Basidiomycota were found to constitute the main members of the fungal community.
Cluster analysis showed a clear separation of different instar larval gut microbiotas, with the exception of 2nd instar larvae, which clustered with early or late instar larval groups because 2nd instar larvae were in a developmental transition from 1st instar to 3rd instar larvae. Additionally, there was also a significant separation of the microbiota community composition in the same instar larval stage collected from different locations. The microbiota of the insect gut is influenced by several factors, such as the environmental habitat, diet, immune system, physiology, and phylogeny of host (Yun et al., 2014;Yao et al., 2016;Martinson et al., 2017;Mistry et al., 2017). Considering that larvae of both populations were developing in different hosts and locations, it is highly likely that the diet and environment played a role in shaping the microbial community structure of larvae. Although the structure of microbiota was different within the two B. minax field populations, a similar shift pattern across the life stages of larvae was identified. This corresponding dynamic change in the gut microbiota during the host development indicates that host physiology may play a vital role in determining the microbiota composition.
Across the different larval developmental stages, we found that only a few bacteria stably resided in the gut of 3rd instar larvae. Most bacterial community members belonged to Enterobacteriaceae, Acetobacteraceae, Lactobacillaceae, and Leuconostocaceae families. Enterobacteriaceae is also the dominant commensal bacterial family in the gut and reproductive organ of B. minax adults (Wang et al., 2014) and in the gut of several other tephritids, including Anastrepha, Ceratitis, and Dacus (Drew and Lloyd, 1991;Kuzina et al., 2001;Behar et al., 2008). Enterobacteriaceae bacterial family have been reported to participate in nitrogen fixation and pectinolysis (Behar et al., 2005(Behar et al., , 2008. Recently, Zhao et al. (2018) found that Enterobacteriaceae, Acetobacteraceae, and Lactobacillaceae were the main families in B. dorsalis larvae and likely participated in sugar metabolism. Previously study showed that bacteria in Acetobacteraceae and Lactobacillaceae families played a vital role in Drosophila larval development by modulating sugar metabolism related insulin pathways (Shin et al., 2011;Storelli et al., 2011). Leuconostoc, Fructobacillus, Oenococcus, and Weissella are representative genera in Leuconostocaceae family. It has been reported that Leuconostoc species can ferment fructose (Ljungdahl, 1962). A specific group of bacteria that solely inhabit fructose-rich niches has been characterized and defined as fructophilic lactic acid bacteria (FLAB) (Endo and Salminen, 2013). FLAB have been found in the guts of bees, tropical fruit flies and giant ants, which live by feeding on flowers or fructose diet (Thaochan et al., 2010;He et al., 2011;Koch and Schmid-Hempel, 2011). It is well known that the sugar content increases as a fruit matures (Chen et al., 2012). Following the oviposition of B. minax, a ripening process begins in the fruit, resulting in sugar-rich fruit (Fischer et al., 2016). This high sugar content provides a favorite substrate for these bacteria to colonize and can be fermented by larvae-associated gut microbiota. In our study, the dynamic increases of three important bacterial families Acetobacteraceae, Lactobacillaceae, and Leuconostocaceae were consistent with the PICRUSt prediction of increased carbohydrate metabolism pathways. Spearman correlation analysis further showed that these three bacterial families displayed positive correlations with carbohydrate metabolism pathways, suggesting that these microbial families can help the host to remarkably adapt to this environment, digest fruit, and provide amino acids or other nutrients. Of note, the relative abundances of these three bacterial families were different in the same instar larval stage developed in different plant hosts, higher in DJK than in ZG samples (Supplementary Table S1). This difference in bacterial abundance may explain why some metabolic pathways, such as fructose and mannose metabolism, were more active in DJK larvae than in ZG. It is plausible that the difference in bacterial abundance was influenced by the presence of different metabolic components in orange and mandarin orange, such as the sugar-acid ratio in fruit. In future, it would be interesting to uncover the relationship among the insect, its gut microbiota and its plant host. There was also a sharp decrease in the relative abundances of several bacterial families during the development of larvae, including Brucellaceae, Chitinophagaceae, Moraxellaceae, and Bacillaceae. These bacterial families may represent transient microbes or may perform their function at a special development stage of larvae. Oxalobacteraceae, Comamonadaceae, Caulobacteraceae, Enterococcaceae, Bradyrhizobiaceae, and Sphingomonadaceae were only detected in the gut of 1st or 2nd instar larvae from ZG, and they decreased to an undetectable level in 3rd instar larvae; thus, it is plausible that environmental microbes that originate from the egg and native plant diet will be specially assembled in newly FIGURE 6 | Comparison of predicted KEGG pathways of larvae and Spearman correlation analysis. The inferred metabolic pathways are shown with carbohydrate metabolism (A), fructose and mannose metabolism (B), pentose phosphate pathway (C), and glycolysis/gluconeogenesis (D). Spearman correlation analysis between the carbohydrate metabolic pathways and bacteria (taxonomic breakdown at the family level) that vary significantly during the development of larvae is presented in heatmaps (E). Multiple comparisons were performed with one-way analysis of variance and Tukey's post hoc test. Different letters indicate a significant difference between different samples (p < 0.05). * p < 0.05 and * * p < 0.01, calculated by Student's t-test. In heatmaps, blue represents negative correlation between gene expression and bacterial genus, and red represents a positive correlation, * P < 0.05. emerged larvae . This large difference of early instar larvae microbiota may explain why the dynamic change patterns of most metabolic pathways across the development stages of larvae were different between the two field populations.
Most research to elucidate insect-microbiota interactions has solely focused on bacteria. However, fungus, as an important member of this interaction relationship, has largely been overlooked (Broderick and Lemaitre, 2012). In the present study, Ascomycota and Basidiomycota were dominant in the mycobiota of B. minax larvae. Several orders of insects are intimately associated with yeasts, including flies, bees, beetles and ants (Becher et al., 2012;Vogel et al., 2017), and most of these associated yeasts belong to the Ascomycota phylum (Chandler et al., 2012). The functions of yeast have largely been studied in Drosophila. Yeasts are vital for larval development because they provide the host with vitamins, sterols and amino acids (Becher et al., 2012;Broderick and Lemaitre, 2012;Yamada et al., 2015) and can attract Drosophila via ester production (Christiaens et al., 2014), inducing Drosophila egg-laying behavior (Becher et al., 2012, Fischer et al., 2016. More recently, Fischer et al. (2016) demonstrated that Drosophila prefers a Saccharomyces-Acetobacter co-culture to the same microorganisms grown separately, and Acetobacter metabolism of ethanol derived from Saccharomyces fungus was essential for this co-culture preference, indicating that the coexistence of this bacterial, and fungal microbiota relationship is beneficial for the host.
In our study, we found that there was a sharp increase in Pichiaceae yeast in accordance with Acetobacteraceae bacteria. The functions of Pichiaceae yeast and whether there is a bacteriafungus interaction within the microbiome of B. minax larvae, as has been shown in Drosophila, remain unknown. Strikingly, there was also a sharp decrease in the relative abundances of several fungi, including Thermoascaceae, Plectosphaerellaceae, Aspergillaceae, Metschnikowiaceae, Lasiosphaeriaceae, Cordycipitaceae, Trichocomaceae, Nectriaceae, Chaetomiaceae, Mortierellaceae, and Russulaceae. These fungi may represent transient microbes or opportunistic pathogenic fungi of citrus. B. minax larvae develop within rotting fruit that may contain these environment fungi, while host metabolites, such as acetic acid, or gut-associated microbes may inhibit these heterotrophic microorganisms during larval development (Kang et al., 2003).

CONCLUSION
In conclusion, the present study demonstrated that the gut microbiota of B. minax larvae greatly differ across different development stages and that several members of microbiota exhibit similar shift patterns in both field populations, which may be indispensable to their adaption to the environment and development of the host. More research is required to explore host-microbe interactions in B. minax. These findings will advance the understanding of the gut microbiota of this harmful citrus pest and may provide clues to develop novel biocontrol strategies against B. minax.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the PRJNA542035.

AUTHOR CONTRIBUTIONS
HZ and ZY conceived and designed the project. ZY performed the experiments and analyzed the data. ZY and ZC generated the graphs. QM, SB, YW, PZ, and HM collected the samples and extracted the DNA. HZ, ZY, QM, and MR wrote the manuscript.  In every panel, the left side represents the population collected from ZG (n = 4), and the right side represents the population collected from DJK (n = 4). Multiple comparisons were performed in the two group separately. Different letters indicate a significant difference between different instars in each group (p < 0.05, one-way ANOVA, Tukey post hoc test). The lines show median values per region window, and the shaded area denotes the estimated 95% confidence interval. Boxplots showing the mean abundance of fungi in the larval gut. In every panel, the left side represents the population collected from ZG (n = 4), and the right side represents the population collected from DJK (n = 4). Multiple comparisons were performed in the two group separately. Different letters indicate a significant difference between different instars in each group (p < 0.05, one-way ANOVA, Tukey post hoc test). The lines show median values per region window, and the shaded area denotes the estimated 95% confidence interval.
FIGURE S6 | Predicted gut bacterial functions of larvae from two field populations. (A-J) the inferred metabolic pathways amino acid metabolism, energy metabolism, metabolism of cofactors and vitamins, xenobiotics biodegradation and metabolism, lipid metabolism, nucleotide metabolism, glycan biosynthesis and metabolism, enzyme families, metabolism of other amino acids and metabolism of terpenoids and polyketides are shown at the second hierarchical level (n = 4). Multiple comparisons were performed in the two group separately. Different letters indicate a significant difference between different instars in each group (p < 0.05, one-way ANOVA, Tukey post hoc test).
TABLE S1 | Comparisons of the relative abundance of gut bacteria from larval and adult samples at family level. DATA SHEET S1 | The shared OTUs across the lifespan of larvae.
DATA SHEET S2 | The relative abundance of KEGG pathways in larvae predicted by PICRUSt.