Composition and function of the microbiotas in the different parts of the midgut of Pyrrhocoris sibiricus (Hemiptera: Pyrrhocoridae) revealed using high-throughput sequencing of 16S rRNA

In pyrrhocorids, digestion of food occurs mainly in the midgut, which is divided into four parts (M1–M4), and takes between three and four days. Food is retained in M1 for about 5 h and passes quickly through M4. However, food is retained in M2 and M3 much longer, about 70 to 90 h. The different stages in digestion may be infl uenced by different microbial populations in the different parts of the midgut. In the present study, the microbiota in the four parts of the midgut of Pyrrhocoris sibiricus were analysed in detail using high-throughput sequencing of the 16S rRNA V3–V4 region. The most abundant bacteria in M3 were Actinobacteria (Coriobacteriaceae) whereas it was Proteobacteria (gammaproteobacteria) in M1, M2 and M4. Actinobacteria was the second most abundant bacterial group in M2. According to the PCA analysis, M2 and M3 have the most similar bacterial communities. Burkholderia, which is closely related to the plant-associated beneficial and environmental (PBE) group, was also found in M1, M2 and M4. P redictive functional profi les of the metagenomes revealed that metabolism mostly occurred in M2 and M3. The PICRUSt results were consistent with the 16S rRNA metagenomic analysis and indicate that the bacteria in M2 and M3 play an important role in degrading complex dietary components. * The fi rst two authors contributed equally to the present study. ** Corresponding author; e-mail: zh_hufang@sohu.com INTRODUCTION Symbiotic interactions between animals and benefi cial microorganisms are very common in nature, including insects (Buchner, 1965; Smith, 1989; Moran et al., 2008). In many cases, these symbionts are housed in the gut of their hosts and play an essential nutritional role in maintaining the fi tness of the host insect. Their role can involve detoxifying plant allelochemicals, promoting digestion and providing digestive enzymes or nutrients, especially when the insect is either a blood, sap or cellulose feeder (Buchner, 1965; Douglas, 1992; Akman et al., 2002; Genta et al., 2006; Sudakaran et al., 2012; Taylor et al., 2014; Onchuru et al., 2018). The Pyrrhocoridae are terrestrial bugs of which there are ~340 species belonging to 33 genera worldwide (Schaefer & Ahmad, 2000; Henry, 2009). Most of them prefer to feed on the seeds of Malvales plants, which are vitamin limited and avoided by other insects due to their phytochemical defences (Allen et al., 1967; Abou-Donia, 1976; Ahmad & Schaefer, 1987; Kristenová et al., 2011). The midgut of PyrEur. J. Entomol. 117: 352–371, 2020 doi: 10.14411/eje.2020.040


INTRODUCTION
Symbiotic interactions between animals and benefi cial microorganisms are very common in nature, including insects (Buchner, 1965;Smith, 1989;Moran et al., 2008). In many cases, these symbionts are housed in the gut of their hosts and play an essential nutritional role in maintaining the fi tness of the host insect. Their role can involve detoxifying plant allelochemicals, promoting digestion and providing digestive enzymes or nutrients, especially when the insect is either a blood, sap or cellulose feeder (Buchner, 1965;Douglas, 1992;Akman et al., 2002;Genta et al., 2006;Sudakaran et al., 2012;Taylor et al., 2014;Onchuru et al., 2018).
The Pyrrhocoridae are terrestrial bugs of which there are ~340 species belonging to 33 genera worldwide (Schaefer & Ahmad, 2000;Henry, 2009). Most of them prefer to feed on the seeds of Malvales plants, which are vitamin limited and avoided by other insects due to their phytochemical defences (Allen et al., 1967;Abou-Donia, 1976;Ahmad & Schaefer, 1987;Kristenová et al., 2011). The midgut of Pyr-by 2% agarose gel electrophoresis and target fragments were purifi ed using an AxyPrep gel extraction kit (Axygen, Union City, CA, USA ). The purifi ed PCR products were then sequenced on an Illumina Miseq platform using 2 × 300 base pairs (bp) pairedend reads (Personalbio, Shanghai, China).

Sequence analysis
Sequencing data was processed using the Quantitative Insights Into Microbial Ecology (QIIME, v.1.8.0) pipeline as previously described (Caporaso et al., 2010). Resulting sequence with average Phred scores lower than 20, ambiguous bases, mononucleotide repeats longer than 8 bp, or those which had a length shorter than 100 bp were removed. Sequences that passed the fi lter without any mismatches were trimmed and merged using FLASH v.1.2.11 (Magoc & Salzberg, 2011), according to the principle of 98% overlap of 19 bases. The resulting sequences of each sample were clustered into operational taxonomic units (OTUs) based on 97% identity using the UCLUST function in QIIME (Caporaso et al., 2010;Edgar et al., 2010). OTU t axonomic classifi cation was conducted using a BLAST search of the representative sequences set against the Greeng enes database (Release 13.8, http://greengenes.secondgenome.com/) (DeSantis et al., 2006) and using the best hit (Altschul et al., 1997). An OTU table was generated to record the abundance and taxonomy of each OTU in each sample. OTUs with a sequence frequency below 0.001% were not included in the analysis (Bokulich et al., 2013). To minimize the difference of sequencing depth across samples, an averaged, rounded rarefi ed OTU table was generated by averaging 100 evenly resampled OTU subsets under the 90% of the minimum sequencing depth for further analysis.

Bioinformatics analysis
Rarefaction curves were calculated at a 97% similarity level in QIIME. Richness (Chao1 and ACE) and diversity (Simpson and Shannon) indices were calculated for the microbiota using an out table in QIIME. Structural variation in the microbiota in the four parts of the midgut was investigated using Weighted Uni-Frac principal coordinate analysis (PCoA) in QIIME (Lozupo ne & Knight, 2005;Lozupone et al., 2007). A maximum likelihood (ML) tree was constructed based on the Kimura 2-parameter model using MEGA X with 1000 bootstrap replicates (Kumar et al., 2018). The function of the gut microbiota was predicted using PICRUS t v1.1.3 Release and Greengenes as the sequenced reference, according to the online protocol (http://picrust.github. io/picrust/index.html) (DeSantis et al., 2006;McDonald et al., 2012;Langille et al., 2013). Firstly, OTUs which mapped to the GG13.5 database with a 97% similarity were selected by QIIME. The OTUs counts were normalized by dividing by the 16S rRNA gene copy numbers from known bacterial genomes in Integrated Microbial Genomes (IMG); normalized OTU tables were of these plants (Zhang, 1985;Matolin & Štys, 1987). In the present study, 16S rRNA metagenomic analysis was used to characterize the microbes that inhabit the different parts of the midgut of P. sibiricus. Furthermore, a PIC-RUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) analysis was used to determine the functional potentials of the bacteria in the M1, M2, M3 and M4 using the 16S rRNA marker gene.

Sample collection
Adult specimens of P. sibiricus feeding on the seeds of Abutilon theophrasti Medicus were collected between September and December 2018 in Jinzhong, Shanxi Province (37. 68°N, 112.75°E). Prior to dissection, they were soaked in 70% ethyl alcohol for 3 min and then washed 3 times with sterile deionized water in order to remove exogenous contaminants. The guts were then dissected under sterile deionized water using sterilized tweezers and eye scissors under aseptic conditions. Each midgut was separated into four parts, M1, M2, M3 and M4 (Fig. 1), which were treated separately. Since there were three replicates for each part and e ach replicate consisted of a mixture of samples from 20 adults regardless of sex, a total of 60 adults were dissected. Five hundred milligrams of each part were collected and immediately frozen in liquid nitrogen for subsequent DNA extraction. The remains of the material of each of the parts were also frozen in liquid nitrogen and preserved at -80°C.

DNA extraction and sequencing
DNA was extracted using a Mag-Bind Soil DNA extraction kit (O mega, Norcross, GA, USA) according to the m anufacturer's in structions with slight modifi cations to improve the effi ciency of the DNA extraction from the gut samples. DNA concentration was measured using a NanoDrop ND-1000 spectrophotometer (Th ermo Fisher Scientifi c, Waltham, MA, USA) and its quality confi rmed using agarose gel electrophoresis. The V3-V4 hypervariable region of the 16S rRNA gene was then amplifi ed using PCR and the universal primers 338F (5'-ACTCCTACGG-GAGGCAGCA-3') and 806R (5'-GGACTACHVGGGT-WTCTAAT-3') (Mizrahi-Man et al., 2013). The reaction solution consisted of: 5 μl of Q5 reaction buffer (5×), 5 μl of Q5 High-Fidelity GC buffer (5×), 2 μl of dNTPs (2.5 mM), 1 μl each of forward and reverse primer (10 uM), 2 μl of DNA template, 0.25 μl of Q5 High-Fidelity DNA Polymerase (5 U/μl) and 8.75 μl of ddH 2 O. PCR conditions consisted of initial denaturation at 98°C for 5 min, followed by 25 cycles of denaturation at 98°C for 30 s, annealing at 52°C for 30 s and extension at 72°C for 1 min, with a fi nal extension at 72°C for 5 min. PCR products were verifi ed then aligned to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The nearest sequenced taxon index (NSTI) was calculated to quantify the quality of the predictions. Heatmap was generated using hemI (Deng et al., 2014). PCA analysis was conducted in STAMP v.2.1.3 (Parks et al., 2014) and contributions of core genera were generated from PICRUSt results using Excel. The functional differences between two-groups were compared using two-sided Welch's t-test in the Majorbio Cloud (cloud.majorbio.com).

RESULTS
The mic robiota in the M1, M2, M3 and M4 of P. sibiricus were characterized using 16S rRNA high-throughput sequencing. A total of 463,717 reads were generated from all of the independent biological treatments and replicates. The average amplicon length of the 16S rRNA variable V3-V4 region was 420 bp. After quality fi ltering and read merging, 405,095 (87.36%) sequences were obtained (105,842, 96,241, 98,594 and 104,418 sequences for M1, M2, M3 and M4, respectively). After assignment using UCLUST and removing OTUs with a relative abundance less than 0.001%, a total of 1,111, 1,060, 496 and 1,217 OTUs were obtained for M1, M2, M3 and M4, respectively. All taxa with relative abundances > 0.001% are cited in Table S3. Rarefaction curves indicate suffi cient sequenc-  Bacterial composition in terms of genera in the four parts of the midgut. Genera with an abundance lower than 1% are referred to as "Others". ing depth for taxonomic classifi cation ( Fig. 2A). For all samples, estimates of ACE richness indicate that M2 had the richest microbiota and the Shannon index suggests that diversity in M4 was signifi cantly higher than in the other parts of the midgut (Table S1).  4. Phylogenetic analysis based on the sequences of Coriobacterium sp., Gordonibacter sp., Enterococcus sp., Burkholderia sp. and unassigned OTUs in the different parts of the midgut of Pyrrhocoris sibiricus. Bootstrap values more than 70% are displayed on the bacterial phylogenies. Abbreviations of clades: BCC -Burkholderia cepacia complex; PBE -plant-associated benefi cial and environmental group; SBE -stinkbug-associated benefi cial and environmental group.

Diversity of bacteria in gut microbiota
At the phylum level Proteobacteria, Actino bacteria, Firmicutes and Bacter oidetes were present in all samples of the midgut of P. sibi ricus. Proteo bacteria was the most abundant in M1 (56.39%), M2 (51.68%) and M4 (68.31%) followed by Firmicutes in M1 and M4 (31.60% and 15.97%, respectively). In M2, Actinobacteria was the second most abundant, account ing for 35.51% of all bacteria. The M3 region was unique in that Actinobacteria dominated followed by Proteobacteria with abundances of 49.58% and 37.11%, respectively (Fig. S1). Principal component analysis revealed that M2 and M3 were most similar in terms of their microb iota (Fig. 2B).
When OTUs with an abundance lower than 1% were removed, 10, 13, 8 and 18 core genera were identifi ed in M1, M2, M3 and M4, respectively (Fig. 3). All OTU s with abundances greater than 1% are cited in Table S2. For the phylum Proteobacteria, the most dominant bacterial group was the gammaproteobacterial genus Enterobacter, which was present throughout the midgut. An unassig ned OTU is closely related to Klebsiella pneumoniae with a 99.17% similarity (AF130982, Fig. 4). The Betaproteobacterial genus Burkholderia was detected in M1, M2 and M4. Accor ding to the maximum likelihood (ML) tree, which was constructed based on the sequence obtained from nextgeneration sequencing as well as representative sequences of in-and outgroup taxa, the Burkholderia detected in the midgut of P. sibiricus closely matched sequences of the PBE clade (Fig. 4).
Within the phylum Actinobacteria, the genus Corioba cterium was recorded in all samples, mostly with high levels of relative abundance in M2 and M3. According to the ML phylogenetic tree, the two Coriobacterium OTUs were closely related to strains of Corioba cterium glomerans (X79048, FJ554832) that occur in the intestinal tract of P. apterus with 97.66% 16S rRNA gene sequence similarity. A Gordonibacter OTU was identified in P. sibiricus that only had a 92.57% 16S rDNA similarity with its clos-est relative, Gordonibacter pamelaeae (AM886059), and 99.26% with Gordonibacter sp. (JX406494) strain in the gut of P. apterus. This OTU was fi rst clustered with Gordonibacter sp. with an overall mean genetic distance of 0.022 and then with cultured strains of the genus Gordonibacter (Fig. 4).
Members of the phylum Firmicutes also occurred in gut micribiota. One unas signed Firmicutes OTU was closely related to a strain in the gut of P. apterus (JX406495, with 98.03% similarity), one was clustered with a cultured strain of Clostridioides mangenotii (MN646967, with 97.28% similarity) and the other two with Enterococcus faecalis (MN749533, with 98.47% similarity) and E. casselifl avus (MG871248, with 98.08% similarity) (Fig. 4).

Predictive functional profi les
In order to determine differences in the functions of the microbio ta in the different parts of the midgut we carried out a functional analysis using PICRUSt. This revealed that 17 level2 KEGG orthology groups (KOs) differentially occurred in the four parts of the midgut, which is depic ted in a heatmap (Fig. 5A). The mean NSTI values were 0.06, 0.04, 0.02 and 0.08 for M1, M2, M3 and M4, respectively. Signifi cant differences in predicted microbial functions (two-sided Welch's t-test, P < 0.05) were detected between M3 and the other parts of the midgut (Fig. 5A, Table S4 ). PCA analysis showed that M2 and M3 have the most similar predicted microbial functions (Fig. 5B). Contribut ions of core genera to the predicted functions, which differ signifi cantly among groups, are presented in Fig.  6. Enterobacter and Burkholderia contribute most to the predicted functions in M1, M2 and M4, while in M3 it is Enterobacter and Coriobacterium (Fig. 6).

DISCUSSION
In the present study, the microbiota inhabiting the midgut of P. sibiri cus were characterized using 16S rRNA metagenomic analysis. The results indicate that it differs in the different parts of the midgut (Fig. S1). Proteobacteria are the dominant bacteria in M1, M2 and M4. In M3, however, there is a distinct microbiota, consisting predominantly of Actinobacteria (Fig. S1). These results are in accord with previous studies on D. fasciatus and P. apterus (Sudakaran et al., 2012;Salem et al., 2013). The compo sition of the microbiota in M2 is particularly interesting in having a similar microbial profile of proteobacteria to M1 whilst sharing some core bacterial taxa with M3 (Fig. 2).
There were, however, some differences, the most notable being the distribution of Gordonibacter sp. and C. glomerans. In previou s studies, they were recorded mainly in M3 (Sudakaran et al., 2012). However, they also make up a high proportion of the microbiota in M2. Furthermore, Lactococcus lactis and an undescribed Rickettsiales bacterium, which are core bacteria in M3 in P. apterus (Sudakaran et al., 2012), were absent in all our samples. We also recorded a relatively high abundance of Burkholderia and Enterobacter. Previously, Burkholderia was thought to be absent in Pyrrhocoridae, including Pyrrhocoris (Sudakaran et al., 2015), however, Gordon et al. (2016) report a high concentration of Burkholderia in the pyrrhocorid Euscopus rufi pes. Hence, particular symbionts may be restricted to certain species of Pyrrhocoridae.
Species of Burkholder ia are symbio tic partners of species of the superfamilies Lygaeoidea and Coreoidea, and particularly the Largidae i n Pyrrhocoro idea (Kikuchi et al., 2011a;Itoh et al., 2014;Gordon et al., 2016). Generally , Burkholderia are acquired directly from the soil by early instar nymphs and stored in well-developed crypts in the posterior midgut (Kikuchi et al., 2011b;Boucias et al., 2012). However, in P. sibiricus, these symbionts are not only to the posterior midgut but also present in M1, M2 and M4. Early microscopic studies indicate that Burkholderia in Hemiptera belong to three distinct clades: the Burkholderia cepacia complex (BCC), the plant-associated beneficial and environmental (PBE) group, and the stinkbug-associated beneficial and environmental group (SBE) (Itoh et al., 2014). Unlike other Burkholderia associations in Alydidae (SBE clade) and Blissidae (all clades), the Burkholderia in P. sibiricus are closely rel ated to the PBE clade (Itoh et al., 2014;Takeshita et al., 2015;Gordon et al., 2016). The function of Burkholderia symbionts i n stinkbugs is unknown (Kaltenpoth & Flórez, 2020). It is hypothesized (Takeshita et al., 2015) that they fi x nitrogen for their insect hosts, however, this is still be tested experimentally. In this study, Burkholderia contributes a high proportion of th e genes in the metagenome associated with xenobiotic biodegradation and m etabolism than any of the other pathways studied. This may indicate that Burkholderia in P. sibiricus provides a secondary context-dependent benefi t, namely resistance to insecticides, which is reported in the bean bug, Riptortus pedestris and oriental chinch bug, Cavelerius saccharivorus (Kikuchi et al., 2012). Many herbivorous insects have symbionts that supply nitrogen, essential amino acids, B vitamins and sterols, which are not present in the plant material they consume (Jones, 1984;Douglas, 1992Douglas, , 1998. Seeds of Malvales are vitamin-limited and are generally avoided by other phytophagous bugs. Furthermore, phytochemical defences (gossypol and cyclopropeno ic fatty acids) of the seeds interfere with digestion and so retarded growth, which can result in sterility (Allen et al., 1967;Abou-Donia, 1976;Kristenová et al., 2011). In previous studies, symbionts of Coriobacteriaceae are reported to activate food polyphenols and dietary phytoestrogens (Clavel et al., 2014;Chew et al., 2018). Some species, such as C. glomerans, can supply B vitamins to the host and protect them from parasites (Salem et al., 2014;Onchuru et al., 2018) and species of Clostridia play a role in the fermentation of carbohydrates, including the degradation of cellulose (Makonde et al., 2013;Sabree & Moran, 2014). Enterobacter belonging to the Enterobacteriaceae are present in the gut of Mediterranean fruit fl ies (Ceratitis capitate) (Aharon et al., 2013;Kyritsis et al., 2019). They are important in nitrogen and carbon metabolism and used as a dietary supplement (probiotic) in the diets to rear their larvae (Behar et al., 2005(Behar et al., , 2008Ben-Yosef et al., 2008;Augustinos et al., 2015;Kyritsis et al., 2019). Our study revealed that Enterobacter, Coriobacterium and Burkholderia make up a big proportion of the bacteria and contribute to most of the functions. The heatmap of the metabolic processes occurring in the different parts of the midgut reveals that M2 and M3 are involved in carbohydrate, amino acid, energy, lipid and vitamin metabolism, which accords with our 16S rRNA metagenomic results, which also indicate that these parts of the midgut are the main sites for food digestion, but the mechanism needs further study.
With the help of bacterial symbionts, Hemiptera are able to exploit plant tissues, such as xylem and phloem, which are nutrient-limited and contain phytochemicals (Buchner, 1965). In recent de cades, our knowledge of Hemiptera and their microbes has increased considerably, particularly the identifi cation of their core microbiota, genome sequencing, phylogenetic relationships, and the role of symbionts in the provision of nutrients, insecticide-resistance and defence against parasites (Kikuchi et al., 2012;Stackebrandt et al., 2013;Kaiwa et al., 2014;Hosokawa et al., 2016;Sudakaran et al., 2017;Onchuru et al., 2018;Onchuru & Kaltenpoth, 2019;Kaltenpoth & Flórez, 2020). However, few studies have focused on the role symbionts in the digestion of food. Although further study is needed to verify the exact mechanism, our results indicate that the microbiota in the M2 and M3 parts of the midgut play an important role in degrading complex dietary components. Supplementary material follows (Tables S1-S4, Fig. S1).