The Impact of Environmental Heterogeneity and Life Stage on the Hindgut Microbiota of Holotrichia parallela Larvae (Coleoptera: Scarabaeidae)

Gut microbiota has diverse ecological and evolutionary effects on its hosts. However, the ways in which it responds to environmental heterogeneity and host physiology remain poorly understood. To this end, we surveyed intestinal microbiota of Holotrichia parallela larvae at different instars and from different geographic regions. Bacterial 16S rRNA gene clone libraries were constructed and clones were subsequently screened by DGGE and sequenced. Firmicutes and Proteobacteria were the major phyla, and bacteria belonging to Ruminococcaceae, Lachnospiraceae, Enterobacteriaceae, Desulfovibrionaceae and Rhodocyclaceae families were commonly found in all natural populations. However, bacterial diversity (Chao1 and Shannon indices) and community structure varied across host populations, and the observed variation can be explained by soil pH, organic carbon and total nitrogen, and the climate factors (e.g., mean annual temperature) of the locations where the populations were sampled. Furthermore, increases in the species richness and diversity of gut microbiota were observed during larval growth. Bacteroidetes comprised the dominant group in the first instar; however, Firmicutes composed the majority of the hindgut microbiota during the second and third instars. Our results suggest that the gut's bacterial community changes in response to environmental heterogeneity and host's physiology, possibly to meet the host's ecological needs or physiological demands.


Introduction
It is well known that a wide diversity of insects harbor a number of commensal and mutualistic micro-organisms that have important ecological and evolutionary effects on their hosts [1,2,3,4]. Previous studies have shown that these microbial communities can provide their hosts with many rewards, including a supply of nutrients [1,4,5,6,7,8], contributions to host reproduction and survival [9,10], protection against natural predators and parasites [11,12], detoxification [13,14], and enhanced social interactions [15]. Although the abundance and distribution of these microbial communities have been recently explored in diverse insects, how environmental factors and host physiology affect the composition and diversity of the intestinal microbiota in the host is mostly unknown [10,16,17]. It has been shown that physicochemical parameters can control the vertical distribution of marine bacterial communities [18,19]. Furthermore, climate change [20], soil attributes [21], and plant species [22] are also important in determining the structure and composition of soil bacterial communities. However, whether and how the environmental conditions influence the structure of the bacterial communities in the guts of soil-dwelling insects remains poorly understood.
Scarab larvae live in soil and feed on plant roots and other organic matter of low nutritive value [23]. Like the bulbous paunch of termites, the hindgut of scarab larvae expands to a structure called the fermentation chamber, which houses dense microbial communities [24,25,26,27,28]. Additionally, previous studies have revealed that large and morphologically diverse microbial communities, as well as numerous novel lineages, are typically found in the intestinal tract of scarab larvae [27,29]. These characteristics suggest that scarab larvae could serve as a useful model to study the relationships between the gut's microbial community and the host. However, scarab larvae are widely distributed around the globe, and the larval feeding ability differs significantly between the first instar and the second and third instars. It is believed that the scarab's dense gut bacteria may contribute to digestion, nutrition and methanogenesis [24,30,31], but questions remain about how ecological factors and the host's developmental stage affect the bacterial community in the gut of the scarab larvae.
The large black chafer Holotrichia parallela is an important peanut pest in China, and it causes great economic losses [32]. The larvae of H. parallela live in the soil and prefer to feed on peanuts. They also have an enlarged bulbous hindgut like other scarab beetles. In the present study, the hindgut bacterial community in H. parallela larvae from different geographic locations and instars were investigated using 16S rRNA clone library construction in combination with Denaturing Gradient Gel Electrophoresis (DGGE). By surveying the bacterial composition and variation of these communities, we examined the following subjects: (i) the composition and diversity of gut bacterial communities of H. parallela larvae; (ii) whether the community composition differs between natural population or larval stages; and (iii) if variation exists, whether and how environmental factors or life stages account for this variation. This information is important for the understanding of the symbiotic relationship between the bacteria and the scarab, as well as the mechanisms that determine gut microbiota composition and dynamics.

Ethics Statement
H. parallela has not been notified under any act or laws and rules thereof of the Government of China or any of the Province governments of Liaoning, Tianjin, Shandong, Hubei, Zhejiang, Fujian, Guangdong, Guangxi, Sichuan, Ningxia as an endangered or threatened species restricting or regulating its collection and observation. No permits were required, for collecting the larvae from the field since all the land accessed is not privately-owned or protected in any way, and H. parallela is not an endangered species affecting the biodiversity status.

Insect rearing and collection
In this study, all animals were handled in strict accordance with good animal practice as defined by the relevant national and/or local animal welfare bodies. A laboratory population of H. parallela was reared in a terrarium filled with organic soil held at 2761uC and on a 14:10 h light-dark photoperiod. Different instar larvae were fed with potatoes on the soil surface. Before dissection, individual larvae were selected and maintained separately in soilfree sterile containers and kept from food for 3 days, to reduce food residue in the gut; however, absorbent cotton soaked with sterile distilled water was placed into the containers and replaced every day.
Third-instar larvae of H. parallela were field-collected from ten different locations throughout China (Table 1, Fig. 1) from July through September 2008. These sampling regions were selected due to their different ecoclimatic characteristics and because the white grubs of H. parallela cause significant damage to the peanut production of these areas. From each region, 20-30 larvae were caught. After returning to the laboratory from the collection sites, all of the larvae were maintained separately in soil-free sterile containers and dissected immediately. Soil samples from the immediate surroundings of the larvae at each site were also collected for analysis.

Insect dissection
Individual healthy larvae were dissected according to the modified method described by Zhang and Jackson [28]. The intestinal tract was carefully removed from the body; the hindgut (average weight 0.16 g gut 21 ; five determinations), from the pyloric sphincter to the rectum including the modified fermentation sac, was placed into a 1.5-ml Eppendorf tube containing 0.5 ml of phosphate buffered saline. The hindgut sections of six individual larvae from each site or instar were pooled before DNA extraction.

Extraction and purification of DNA
After mixing well, the hindgut samples were homogenized and centrifuged (130006g for 5 min), and then the hindgut pellet was used for the DNA extraction. The DNA extractions were carried out using the enzymatic lysis protocol described in Yang et al. [33] with the following modifications: each pellet was suspended in 557 ml of Tris-EDTA (TE) buffer (100 mM Tris-HCl pH 7.5, 50 mM EDTA) and 10 ml lysozyme (238 mg mL 21 ) was added. After incubating at 37uC for 20 min, 3 ml proteinase K (20 mg mL 21 ) and 30 ml 10% sodium dodecyl sulfate were added. The mixture was then incubated for 1 h at 37uC. After the incubation step, 100 ml NaCl (5 mol L 21 ) and 80 ml CTAB/NaCl were added, followed by an incubation at 65uC for 10 min. The supernatant was then extracted using consecutive phenol-chloroform-isoamyl alcohol (25:24:1 by volume), isopropanol, and ethanol precipitation steps. Finally, the DNA was re-suspended in 100 ml of TE buffer and stored at 220uC until it was used for Polymerase Chain Reaction (PCR) and further analysis.

PCR amplification of 16S rRNA gene
The variable V6-V8 region of the 16S rRNA gene was amplified using the DGGE primers F-968-GC and R-1401 [34]. The PCR was performed as previously described in Penders et al. [35], except that the annealing temperature was 58uC. Since PCR-based DNA techniques may introduce bias during PCR amplifications, it's very important to minimize PCR artifacts in clone library construction when study an environmental sample with high complex microbial diversity [36,37]. In fact, we had performed a preliminary experiment to test whether a ''reconditioning PCR'' step can reduce the presence of heteroduplex sequences in our gut samples. After analyzing by DGGE, a weak band with novel DGGE mobility was observed in the products of common PCR protocol, But in products of PCR with reconditioning PCR step, no novel visible band was observed. So, ''reconditioning PCR'', was performed as suggested by Acinas et al. [37] and Thompson et al. [38], to reduce PCR-produced biases and artifacts. After amplification, PCR products were checked by

Construction of 16S rRNA gene libraries
Triplicate PCR products were pooled and cleaned using an EasyPure Quick Gel Extraction Kit (BeiJing TransGen Biotech Co., Ltd., BeiJing, China) following the manufacturer's instructions. Amplified gene fragments were cloned using the pEASY-T1 clone kit (BeiJing TransGen Biotech Co., Ltd., China) with blue-white selection. Approximately 90 white colonies were obtained from each of the three larval stages, and 170 white colonies were chosen from each natural population. Plasmid DNA from the cultured transformants was isolated by alkaline lysis [39], followed by PCR amplification using the M13 forward and reverse primers provided by the kit following manufacturer's instructions (pEASY-T1 clone kit, BeiJing TransGen Biotech Co., Ltd., China). For each clone library, positive clones that containing inserts of approximately 430 bp were selected and subsequently screened by DGGE.

Screening of the 16S rRNA gene clone library by PCR and DGGE
The DGGE approach is widely used for monitoring bacterial growth and analyzing bacterial communities [28,40,41]. However, limitations of DGGE have also been pointed out by various authors, such as only bacterial populations that make up more than 1% of the total community can be detected [42,43,44], a complex community that comprised of numerous populations in relatively equivalent proportions will result in a smear of bands, which makes it difficult to identify individual populations and excise DGGE bands [45]. In order to overcome these deficiencies, we constructed the 16S rDNA clone libraries and screened the different clones by DGGE. This strategy combines the advantage of the original cloning method of the 16S rRNA gene with the advantage of DGGE, and allows screening different clones and analyzing the structure of the bacterial community in one gel.
DGGE was performed using the DCode mutation detect system (Bio-Rad, USA). One-microliter PCR product of each positive clone that amplified by the M13 primer set in the above step was served as template DNA in a second PCR step using the F-968-GC and R1401 primers as described previously. Then 20 ml samples of the obtained PCR products were subsequently loaded on to 8% (w/v) polyacrylamide gels with a denaturing gradient of 30-70% (100% denaturant solution consists of 7 M urea and 40% formamide) and run for 20 h at 80 V at a constant temperature of 60uC in 0.56 Tris-Acetate-EDTA buffer (pH 7.4). To standardize the DGGE gels, a mixture of nine cloned DNA fragments with different electrophoretic positions was loaded as a marker. After electrophoresis, the DGGE bands were visualized with silver nitrate [46], and the DGGE gels were analyzed by GELCOM-PARII software(Applied Maths NV, Belgium). Each clone library was analyzed individually, and clones with the same electrophoretic mobility in the gel were temporarily grouped into one group (Fig. 2). Once DGGE analysis was finished, the number of different groups in each clone library was recorded. The number of clones in each group was also recorded.

Sequencing and Sequence analysis
For each clone library, a representative clone from each DGGE group was sequenced. All sequences were checked for chimeric artifacts using the check-chimera program of the Ribosomal Database Project (RDP) [47], and aligned against those found in the GenBank nr database, in the RDP II database [47], and on the EzTaxon server [48] using the Basic Local Alignment and Search Tool (BLAST) algorithm [49]. These nonchimeric sequences data have been submitted to the GenBank database under accession numbers JF964265-JF964858 and JN006162-JN006258.
All the sequences from the bacterial clone libraries were aligned with ClustalX(http://www.clustal.org/),and distance matrices were calculated according to the Jukes-Cantor distance Model in the PHYLIP 3.69 software package (http://evolution.gs. washington.edu/phylip.html) to group the sequences into operational taxonomic units (OTUs). Sequences were clustered by the program MOTHUR into OTUs using the furthest neighbor method and identity cut-offs of 0.03 representing approximately species-level clustering [50]. The OTU was also represented by the name of the representative clone in further phylogenetic analyses. In each clone library, if several groups were classified as the same OTU finally, the number of clones in these groups were added up to calculate the number of clones that this OTU contained. The

Phylogenetic analysis
Phylogenetic analysis of these partial 16S rRNA gene sequences was performed using the MEGA4.0 software [51]. For phylogenetic analysis, a representative clone of each OTU was chosen from the clone libraries, and their sequences were aligned with the nearest neighbors found in the NCBI database and SILVA web server [http://www.arbsilva. de/] [52]. Clones obtained from the larval gut of other scarab beetles that were available in NCBI database were also added to the alignment. Neighbour-joining (NJ) analysis was performed by using MEGA4.0 [51]. The NJ tree was constructed from the distance matrix calculated by the algorithm of Kimura's two-parameter model. Bootstrap confidence values were obtained with 1000 resamplings.

Soil analyses and climatic data
Air-dried soil samples that had been sieved using a 2.0-mm sieve were used for analysis. Soil pH values were determined using a water (water: soil = 2.5:1) suspension [53]. The total amount of organic carbon in the soil (hereafter ''Total soil C'') was analyzed using the H 2 SO 4 -K 2 Cr 2 O 7 oxidation method with an Alpkem autoanalyzer (Kjektec System 1026 Distilling. Unit, Sweden). The total amount of nitrogen in the soil (hereafter ''Total soil N'') was measured using the Kjeldahl acid-digestion method. All of these analyses were completed in the Key Laboratory of Subtropical Agriculture and Environment, Ministry of Agriculture, in the Huazhong Agriculture University. The mean annual temperature (Annu), January low temperature (LowT), monthly temperature range (Ave, defined as the 12-month average of the differences between the monthly mean maximum and minimum temperatures), and mean annual precipitation (PRE) climatic data represent 40-year averages recorded between 1961 and 2000 at various locations throughout China. These data were obtained from http://cdc.cma.gov.cn, which is published by the China Meteorological Administration.

Statistical analyses
The richness estimations and diversity indices were calculated using EstimateS (Version 8.2, http://viceroy.eeb.uconn.edu/ EstimateS) software. The bias-corrected Chao1 estimator of species richness was calculated after 1,000 randomizations of sampling without replacement. The diversity of the sampled sequence set was estimated using the Simpson and Shannon indices (H9) in the EstimateS application. The coverage of the clone library is given as: C = 12(n1/N), where n1 is the number of clones that occurred only once in the library, and N is the total number of clones examined [54]. Rarefaction curves were produced using the Analytic Rarefaction 1.3 program, which is available online at http://www.biology.ualberta.ca/jbrzusto/ rarefact.php.
The correlation between the Chao1 estimates and environmental factors was analyzed using the Spearman nonparametric correlation test available in the SPSS 16.0 software (SPSS Inc., Chicago, IL).
#-LIBSHUFF, a tool in the mothur software package that implements the Cramer-von Mises test statistic, was used to compare the clone libraries and determine the degree of similarity between them [50]. The relationship between the hindgut microbial community composition and environmental factors was analyzed by redundancy discriminate analysis (RDA) using CANOCO software (Microcomputer Power, Ithaca, New York). The data matrices containing the species data were log(x+1) transformed before analysis. Out of all of the environmental variables, the environmental factors that best described the most influential gradients were identified by manual forward selection and subjected to further analysis. We used a Monte Carlo permutation test based on 499 random permutations to test the significance (P,0.05) of the relationship between the explanatory variables and the community composition.

Phylogenetic analyses of the dominant hindgut bacteria of H. parallela
A total of ten clone libraries were constructed to represent the bacterial communities naturally present in the hindgut of scarab beetle larvae. After DGGE analysis, 44-78 representative clones from each clone library (total 594 clones) were selected for sequencing ( Table 1). None of the clones was discarded as chimeric, and 408 sequences (68.7% of the 594 sequenced clones) shared ,97% sequence identity with published 16S rRNA gene sequences in the GenBank database. The 594 representative sequences were identified as Firmicutes, Alpha-, Beta-, Gammaand Deltaproteobacteria, Bacteroidetes, Actinobacteria and Fusobacteria sequences. At a .97% sequence identity threshold, 205 OTUs could be defined. Most of the OTUs belong to Firmicutes (166 OTUs), Deltaproteobacteria (12 OTUs), Gammaproteobacteria (8 OTUs), and Betaproteobacteria (7 OTUs). The remaining groups were represented, to a much smaller extent, by five OTUs for Bacteroidetes, three  OTUs for Actinobacteria, three OTUs for Alphaproteobacteria, and only one representative OTU for Fusobacteria.
Members of the largest group, Firmicutes, mostly belonged to the Clostridia class (Fig. 3), followed by Bacilli and Erysipelotrichi (Fig. 4). Most of the sequences belonging to Clostridia were affiliated with the Lachnospiraceae, Ruminococcaceae, Clostridiaceae, Eubacteriaceae, Veillonellaceae, and Christensenellaceae families (Fig. 3). The largest cluster was Lachnospiraceae, and clones in this cluster were closely related to several clones obtained previously from larvae of scarab beetle Pachnoda ephippiata, Dermolepida albohirtum, and Costelytra zealandica. Also, several clones were affiliated with Clostridium jejuense as the closest relative. Clones in another important cluster (Ruminococcaceae) was closely related to Oscillibacter valericigenes, and also comprised several clones obtained previously from the gut of the larvae of other scarab beetles (Fig. 3). The Clostridiaceae clones formed a distinct lineage among the Clostridia that contained no cultivated representatives. Their closest relatives were clones from other scarab beetles, termite guts and Spodoptera littoralis. Almost all of the OTUs of Deltaproteobacteria were grouped into the Desulfovibrio genus, and their closest relatives were Desulfovibrio cuneatus and some clones from scarab beetle guts (Fig. 5). The OTUs of Gammaproteobacteria were mostly related to the Enterobacteriaceae family, mainly the Klebsiella, Enterobacter and Raoultella genera (Fig. 5). All of the OTUs of the Bacteroidetes were related to the Porphyromonadaceae family. Several clones in this group were affiliated with Dysgonomonas hofstadii as the closest relative, while numerous other clones were closely related to clones from termite guts and other scarab beetles (Fig. 6). One Actinobacteria OTU that belonged to the Promicromonosporaceae family was identified as Promicromonospora pachnodae (99% sequence similarity), a dominant (hemi)cellulolytic bacterium found in the hindgut of the scarab beetle Pachnoda marginata (Cazemier et al., 2003) (Fig. 6).

Diversity analysis and structural comparison of gut microbiota
We evaluated the species richness and diversity of the bacterial communities in the hindgut of larvae of H. parallela from different sampling sites. Coverage values (ranging from 86.1% to 95.6%) and rarefaction curves (Fig. S1) indicated that most prevalent bacterial groups in the hindgut were identified. The number of OTUs identified from the different populations varied from 35 to 49 OTUs (Table 1). Species richness also varied across the host locations, as measured by the Chao1 estimator. The bacterial community from the TJ population (38.8) possessed the least species richness, while the FJ (86.5) population had the highest species richness. In addition, the Shannon diversity index (H9) in each clone library was high (ranging from 3.2 to 3.54), indicating a high diversity of bacterial phylotypes in the hindgut of H. parallela.
One prominent feature of the gut microbiota of H. parallela larvae was the high prevalence of Firmicutes and Proteobacteria sequences (Fig. 7). Sequences affiliated with Firmicutes (mainly Clostridia) were detectable in all of the natural populations, and they were the most prominent group in the hindgut microbiota of the ten populations (ranging from 52% to 86.2%). Sequences related to Beta-, Gammaand Deltaproteobacteria were also present in all natural populations, although at relatively low proportions (Fig. 7). However, even bacteria belonging to these major groups above, such as Ruminococcaceae, Lachnospiraceae, Enterobacteriaceae, and Desulfovibrionaceae, were commonly found in all of the natural populations, and constituted dominant and stable populations in the scarab gut, the composition and complexity of these dominant bacterial groups varied significantly among the ten populations. Only six OTUs (OTUNXn56, OTUNXn53, and OTUGD108 in the class Clostridia, OTUTJV63 in the class Gammaproteobacteria, OTUNXm47 in the class Betaproteobacteria, and OTULN154 in the class Deltaproteobacteria) were shared between the ten natural populations (Fig. 3 and Fig. 5); However, the remaining OTUs were unevenly spread across the host populations, with many appearing in only 1 host population among the ten sampled. In addition, bacteria belonging to other phyla (Bacteroidetes, Actinobacteria, and Fusobacteria) were also unevenly distributed across the host populations (Fig. 7). For example, Bacteroidetes were absent from the GX, HB, and LN populations, whereas bacteria belonging to the Fusobacteria were only detected in the SC population. The similarity indices calculated using the #-LIB-SHUFF tool showed that the communities of the ten clone libraries are significantly different from each other (P,0.0011), with the exception of the clone library of the ZJ population, which showed no significant difference from the FJ and NX populations ( Table S1).

Influence of environmental factors on the hindgut bacterial community
The climate parameters and soil characteristics of the various locations sampled in this study are listed in Table 2. They revealed significant differences among the sampling regions. Positive correlations were found between the Chao1 estimate of richness and the Annu, PRE, and LowT measurements (Methods , Table 3). However, a negative correlation between the Chao1 estimate of richness and the soil pH (r s = 20.739, p = 0.015) was also observed. The Chao1 estimate of richness was not correlated with other environmental factors, including Ave, total soil N, total soil C and the C:N ratio ( Table 3).
The RDA revealed that the bacterial community composition was related to environmental factors, including Ave, Annu, soil pH, total C and total N. The RDA yielded two main axes that explained 22.6% and 18.2%, respectively, of the total variation in bacterial community structure (Fig. 8). Biplot scaling of the RDA further indicated the influence of the different variables on the various bacterial phylotypes. The total C and total N of the soil correlated with the presence of the Ruminococcaceae family, Eubacteriaceae family and Rhodocyclaceae family. Correlations were also detected between the Ave and the Clostridiaceae family, Veillonellaceae family, and family Enterobacteriaceae. Additionally, the soil pH was correlated with the presence of the Raoultella genus, and Annu was correlated with the presence of the Desulfovibrio genus and Klebsiella genus. Additionally, there was a negative correlation of total N with the presence of bacteria from the Lachnospiraceae family, and the presence of the Christensenellaceae family were negatively correlated with Ave.

Gut bacteria community in the larvae of different instars
Besides comparing the bacterial 16S rRNA gene sequence libraries recovered from ten natural populations, we also determined the gut microbial communities in the larvae of different instars. The purposes was to identifying the predominant Overall, three 16S rRNA gene sequence libraries were constructed from all three instars. Approximately 30 representative clones from each clone library were sequenced and analyzed (Table 4). Fifty-eight OTUs over a broad diversity of bacteria were identified, including Firmicutes, Proteobacteria, Actinobateria, and Bacteroidetes ( Fig. 9 and Fig. 10). Twenty-nine OTUs grouped with the Clostridia, and several clones in this group were affiliated with Clostridium nexile, Ruminococcus gauvreauii, and Robinsoniella peoriensis as the closest relative. Other Clostridia clones formed distinct lineages with no cultured representatives, and were closely related to the clones obtained from P. ephippiata and Protaetia brevitarsis (Fig. 9). OTUs affiliated with Proteobacteria were also dominant in the analyzed clones. Ten OTUs were grouped with the Proteobacteria (2 Alphaproteobacteria, 3 Betaproteobacteria, 4 Gammaproteobacteria, and 1 Deltaproteobacteria OTU) (Fig. 10). Nine OTUs were grouped within  the Bacteroidetes (mainly Porphyromonadaceae), and seven OTUs within Bacilli. The remaining three OTUs were grouped within the Actinobacteria (Fig. 10).
Furthermore, significant changes in the complexity of the bacterial population during larval development were also observed ( Fig. 9 and Fig. 10). The first instar larvae harbored a simpler bacterial community than the second and third instar larvae. Bacteroidetes were the mostly commonly found bacteria in the hindgut of the first instar larvae (48 clones, 54.55% of the first instar total clones) (Fig. 11), with one phylotype (represented by clone B77) predominantly found (37 clones, 42.05%). This OTU was also present with a low prevalence in the second-and third instar. Bacterial communities in second and third instar larvae were more complex, and they were dominated by the Ruminococcaceae, Lachnospiraceae, Enterococcaceae and Enterobacteriaceae families, whereas the Bacteroidetes phylum became a minor component of the community (Fig. 11). Overall, four phylotypes, including one Porphyromonadaceae (OTUB77), one Pseudomonadaceae (OTU F35), one Rhodocyclaceae (OTUF2), and one Ruminococcaceae (OTU F26), were present in the hindgut during all three instars, but at different frequencies. However, another 12 phylotypes were also shared by the second and third instar larvae, suggesting the compositions of the communities in the second and third instar larvae were similar to each other.
Although there are several OTUs in common between all three instars, the #-LIBSHUFF analysis revealed clear differences in the microbial community composition between the gut microbiota of first instar larvae and that from the second or third instar larvae Abbreviations are used as follows: Annu, Mean annual temperature; Ave, monthly temperature range; PRE, mean annual precipitation; LowT, January low temperature; Total N, the total amount of nitrogen in the soil; Carbon, the total amount of organic carbon in the soil. doi:10.1371/journal.pone.0057169.t002  (P,0.0083; Table S2). No significant differences were detected between the second and third instars. Chao1 estimator analysis (Table 4) and the rarefaction curves (Fig. S2) both showed that the hindgut bacteria communities in the second and third instar larvae were more diverse than those in the first instar larvae.

Discussion
Using the DGGE technique and the analysis of the 16S rDNA clone libraries, we have identified bacteria that are associated with the H. parallela larval host over a broad geographic distribution. Diverse bacteria belonging to the Firmicutes, Proteobacteria, Bacteroidetes, Actinobacteria and Fusobacteria phyla were identified in the hindgut of H. parallela larvae, with the first two phyla dominating in all of the sampled populations. The microbial distribution of major phyla in the hindgut of H. parallela larvae was similar to that found in other scarabs, such as Melolontha melolontha [26], D.
albohirtum [27], C. zealandica [28], and P. ephippiata [29]. Furthermore, many phylotypes from the hindgut of H. parallela larvae were closely related to the clones obtained from other scarab beetles (P. ephippiata; D. albohirtum, and C. zealandica), and some bacterial species even constituted monophyletic clusters (scarab beetle clusters) together with clones from other scarab beetles. These similarities suggest that although the bacterial community in the intestine of different scarab species is highly diverse, the phylogenetic placement of the taxa from the Clostridia, Actinobacteria, Deltaproteobacteria is correlated with each other. Sometimes these correlations are significant [27], and perhaps the presence of the taxon-specific lineages is indicative of cospeciation between gut microbiota and the host [29].
One prominent feature of the gut microbiota of H. parallela larvae was the high prevalence of Firmicutes sequences, which is similar to previous studies on the gut microbiota of other scarab beetles M. melolontha [26], C. zealandica [28], and P. ephippiata [29]. Firmicutes have been broadly found in the intestine of iguanas and mammals, and these bacteria are thought to have important and core roles in the host's metabolism [55,56]. Microbial symbionts, especially the clostridial bacteria, are responsible for the metabolism of cellulose in the bamboo diet of pandas [57]. This situation also appears to apply to the Scarabaeidae family. Firmicutes contain many cellulolytic Bacillus and Clostridium species. The cellulose and hemicelluloses in roots that are ingested by the scarab larvae can first be hydrolyzed by Firmicutes to their constituent hexoses and pentoses, which can be further fermented to volatile fatty acids (VFA) or other by-products usable by methanogenic bacteria [30]. Additionally, Clostridiaand Bacteroides-related bacteria are responsible for the formation of butyrate, a characteristic product in the hindgut contents of P. ephippiata [31]. Functions have been suggested for some of the other bacterial groups detected in this study. Members of the Enterobacteriaceae family have involved in nitrogen and carbon metabolism in fruit fly Ceratitis capitata, and had an indirect contribution to host fitness [10]. A recent study showed that Desulfovibrio spp. were the most prominent group of microorganisms located almost exclusively in the gut periphery of M. melolontha, and played an important physiological role of sulfate metabolism of the host [26]. Interestingly, in this study, six bacterial species were consistently detected in all of the natural populations, suggesting the existence of a symbiotic relationship between these bacteria and the H. parallela larvae. This evidence further suggests that hindgut bacteria may fulfill important primary functions in the scarab beetle.
The diversity and composition of the hindgut bacterial communities varied significantly between larvae from different geographic location. The bacterial species richness significantly correlated with LowT, Annu, PRE and soil pH. Additionally, the composition and complexity of those dominant bacterial groups were also significantly affected by Ave, Annu and soil characteristics (i.e., pH, total C, total N), suggesting that the environment inhabited by the scarabs can influence the diversity of their gut  bacterial communities. Our results are similar to those found in previous studies [17,58,59,60]. Tsuchida et al. [58] reported a significant positive correlation between the pea aphid U-type symbiont (PAUS) infection and the mean annual temperature and mean annual precipitation. Previous research showed that the prevalence of symbionts in the stinkbugs Acrosternum hilare and Murgantia histrionica was significantly influenced by temperature [60]. Zouache et al. [17] found that habitats can affect the bacterial diversity of wild Aedes populations. Because insects are smallbodied ectotherms, their life cycle duration, survival, and behavior are dependent on ambient temperature [61,62]. A warm climate might favor the survival and reproduction of bacteria in the insect's gut. Therefore, it is appropriate that bacterial species richness was significantly correlated with LowT and Annu. In addition, scarab larvae live in the soil where they feed on organic matter. The soil characteristics and the soil microbial community may also affect the bacterial composition of scarab larvae's gut, and they may have caused the spatial changes in the bacterial communities that were seen in this study. Previous studies on the gut microbes of sand fly Phlebotomus argentipes [63], zebrafish [64] and the symbionts of the stinkbug Riptortus clavatus [65] have shown a strong correlation between insect-microbial associations and the environment in which they reside. Furthermore, it has been shown that soil chemistry characteristics [66,67], soil moisture, and soil temperature [68,69] markedly affect soil bacterial composition and diversity. Previous studies have also shown that plant species and agricultural management regime can influence the activity and composition of soil bacterial communities [22,41,70]. Although most of the sampling sites were peanut fields, the management regime of these fields were different. For example, wheat-corn rotation was the main agricultural management regime where the NX population was collected, but the field that the GX population was collected has been exposed to spring soybean-peanut rotation for several years. Moreover, the GD population and ZJ population were collected from fields subjected to rice-peanut rotation. The difference of agricultural practices in these sampling sites can also have impacts on the bacterial community structure in the soil and were responsible for the bacterial community variations we recorded. However, in this study, we found that some environmental factors that can affect the soil bacterial community, such as soil pH, total C and mean annual precipitation (which can affect soil moisture), also affected the constitution and diversity of the bacterial community of the scarab larvae's gut. This result not only indicates that some of the bacterial species that colonize the gut of scarab beetles might originate from the surrounding soil, but it also highlights, as suggested by Zouache et al. [17], the importance of considering environmental factors (i.e., ecological niches and ecoclimatic characteristics) when analyzing the symbiotic microbiota associated with wild animal populations.
Our results also suggest that the host developmental stage affects the composition of the hindgut bacterial community. The first instar larvae of H. parallela have microbiota with lower diversity, and the hindgut microbiota is dominated by Bacteroidales. However, in the second-and third instar, there were more diverse bacteria species, and the hindgut microbiota during these instars are dominated by several phylotypes such as Clostridiales and Lactobacillales. Such shifts in the dominant gut bacterial species between different life stages have also been detected in Macrotermes gilvus [71], Anopheles gambiae [72], Frankliniella occidentalis [73], and Wohlfahrtia magnifica [74]. The amount of feeding is usually low during the first instar, and it increases during larval development or body size increase [75,76]. The variation in the gut bacterial community might be due, to some extent, to the food source and the digestion process, as suggested by Hongoh et al. [71]. The second and third instar larvae require greater amounts of food than first instar larvae. The increased number of fermentative bacteria groups, like Lactobacillales and Clostridiales, during the second and third instar larvae indicates that these bacteria might be involved in the fermentation process of the food in the gut of the larvae. Alternatively, the animal gut pH, redox conditions, and digestive enzymes present also have significant effects on the gut microbial community [1]. Therefore, it is possible that the hindgut bacterial community shift correlates with changes of gut physiology during larval development. Table S1

Supporting Information
(DOC)