Dynamics of the hindgut microbiota of the Japanese honey bees (Apis cerana japonica) throughout the overwintering period

TL;DR

This study examines how the gut bacteria of Japanese honey bees change during winter and identifies potential new bacteria that may help them survive.

Contribution

The study identifies potential novel bacterial species in the hindgut microbiota of Japanese honey bees and their abundance changes during overwintering.

Findings

The hindgut microbiota of Japanese honey bees is dominated by six core genera, some of which may represent novel species.

Bifidobacterium, Bombilactobacillus, and Lactobacillus were more abundant during overwintering.

The microbiota composition of Japanese honey bees is similar to other Apis cerana populations, suggesting phylogenetic influence.

Abstract

Honey bees play crucial roles as pollinators in natural, agricultural, and ecological systems. The role of gut microbiota in the overwinter survival of honey bees is gaining attention. Compared with Western honey bees (Apis mellifera), Eastern honey bees (Apis cerana) are more tolerant to low-temperature stress. This study compared the hindgut microbiota of the Japanese honey bees (Apis cerana japonica), a subspecies of A. cerana, during the overwintering period (December) with that before overwintering (October) and after overwintering (March) to estimate beneficial hindgut bacteria contributing to survival during the overwintering period. Overall, the hindgut microbiota of A. c. japonica was occupied by Actinobacteriota, Bacteroidota, Firmicutes, and Proteobacteria at the phylum level and Apibacter, Bifidobacterium, Bombilactobacillus, Gilliamella, Lactobacillus, and Snodgrassella at the genus level. The hindgut microbiota composition of A. c. japonica was similar to that of A. cerana in other regions, suggesting that phylogeny influenced the composition. Many sequences assigned to these six core genera showed <98.7% similarity to type strains, indicating potential novel bacterial species. The relative abundance of Bifidobacterium, Bombilactobacillus, and Lactobacillus was higher during overwintering than in other periods. Our findings highlight changes in the core bacteria of the hindgut microbiota of A. c. japonica during overwintering and also suggest the presence of novel candidate bacterial species. The roles of the bacteria that were increased during the overwintering period require further elucidation.

Figures

Hindgut microbiota composition of the Japanese honey bees (Apis cerana japonica) sampled for each period at the (A) phylum and (B) genus levels.BO, before overwintering; OW, during overwintering; AO, after overwintering.

Nonmetric multidimensional scaling (NMDS) ordination plots of hindgut microbiota of the Japanese honey bees (Apis cerana japonica) at three sampling periods.The plot was generated with the Bray–Curtis dissimilarity index based on the ASVs obtained from each sample. BO, before overwintering; OW, during overwintering; and AO, after overwintering.

Full text

Introduction

The gut microbiota of honey bees plays a critical role in their health (Raymann & Moran, 2018), including decomposing dietary compounds (Engel, Martinson & Moran, 2012), producing short-chain fatty acids (SCFAs) as an energy source (Zheng et al., 2017), degrading potentially toxic plant metabolites (Motta et al., 2022), inhibiting the growth of honey bee pathogens (Wu et al., 2014), and stimulating the immune system (Kwong, Mancenido & Moran, 2017; Motta & Moran, 2024). Disruption of the gut microbiota composition due to antibiotic treatment and pesticide exposure causes dysbiosis, leading to host mortality (Raymann, Shaffer & Moran, 2017; Motta & Moran, 2024).

During winter, honey bees survive the severe cold environment in a metabolically and physically active state that is essential for ensuring the colony’s survival until the following spring (Moeller, 1977; Doeke, Frazier & Grozinger, 2015). Cold stress is a major cause of individual and colony mortality in honey bees and also increases the risk of disease and infection outbreaks (Doeke, Frazier & Grozinger, 2015; Xu et al., 2017). During overwintering, feeding is essentially limited to food stored within the colony (pollen, bee bread, and honey). To cope with the surrounding cold stress, the honey bees must maintain the temperature of the colony’s outer edge and core by vigorously vibrating their flight muscles to generate heat (Doeke, Frazier & Grozinger, 2015). The gut microbiota of overwintering honey bees has been increasingly recognized for its beneficial role in survival during the overwintering period. Most studies have focused on the Western honey bees (Apis mellifera), reporting an increase in gut bacteria abundance during this period (Kešnerová et al., 2020) along with changes in gut microbiota composition (Bleau et al., 2020; Kešnerová et al., 2020; Liu et al., 2021; Castelli et al., 2022; Li et al., 2022; Brar et al., 2025). These findings suggest that gut bacteria may play crucial roles in energy absorption and immune function, thereby substantially contributing to survival during the overwintering period.

The Apis genus is naturally distributed across Asia, Europe, and Africa (Ji, 2021). The Western honey bees are widely distributed worldwide, including their native regions of Europe, Africa, and the Middle East, while Eastern honey bees (Apis cerana) are found throughout South, Southeast, and East Asia (Ji, 2021). Compared to Western honey bees, Eastern honey bees exhibit superior cold tolerance and are better adapted to surviving the harsh overwintering period (Li et al., 2012; Xu et al., 2017). In consideration of the advantageous role of the gut microbiota in the overwintering process of honey bees, we hypothesized that the gut microbiota plays a significant role in enabling Eastern honey bees to tolerate cold stress and successfully overwinter. Characterizing the hindgut microbiota of overwintering honey bees will help narrow down the candidate bacteria beneficial for survival under cold and harsh environments, providing novel insights into the symbiotic relationships between honey bees and their hindgut microbiota.

This study focused on the Japanese honey bees (Apis cerana japonica), a subspecies of the Eastern honey bee that is native to Japan. The aim was to develop a comprehensive inventory of the hindgut microbiota of A. c. japonica using high-throughput sequencing of the V3–V4 region of the bacterial 16S rRNA gene. In addition, we compared the composition of the hindgut microbiota before, during, and after the overwintering period to identify distinctive microbial features associated with successful overwintering.

Materials and Methods

Sample collection

The study samples were collected from four A. c. japonica colonies in Ibaraki, Japan. Two colonies were maintained by our laboratory at the National Institute of Environmental Studies in Tsukuba City, while the other two were managed by local beekeepers in Tsukuba City and Inashiki District, respectively. We sampled 30 foragers from each colony using a net with clean plastic cups over three periods: October 2022 (before overwintering, BO group), December 2022 (during overwintering, OW group), and March 2023 (after overwintering, AO group) (Table 1). The total number of samples was 360 (30 honey bees per colony × four colonies × three periods). All samples were immediately placed on ice after collection and stored at −80 °C until DNA extraction.

Table 1: Detailed information on sampling in this study.

DNA extraction

After thawing the honey bees on ice, they were surface-sterilized by immersing them in 70% ethanol for 30 s, followed by rinsing with ultrapure water for 30 s. The hindguts, including the pylorus, ileum, and rectum, were carefully removed on ice using sterile forceps. Ten hindguts were pooled into 2.0 mL sterile tubes containing TE buffer (10 mmol L^−1^ Tris-HCl and 1 mmol L^−1^ EDTA-2Na, pH 8.0) supplemented with 5% (v/v) Triton X-100 (MP Biomedicals, Irvine, CA, USA) and glass beads (1.0 mm diameter). The hindguts were disrupted by three cycles of crushing at 3,200 rpm for 30 s using Beads Crusher µT-12 (Taitec, Saitama, Japan) and 30 s of cooling on ice. The homogenates were centrifuged at 6,000× g for 10 min to sediment debris. Total bacterial DNA was extracted from 180 μL of the resulting supernatant using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacture’s instructions. The DNA concentration of the 36 samples (three replicates per colony × four colonies × three sampling periods) was measured using a NanoDrop One spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

High-throughput sequencing of the V3–V4 region of the bacterial 16S rRNA gene

The bacterial V3–V4 hypervariable regions of the 16S rRNA gene were amplified using a universal primer set comprising 341F (5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCT-NNNNN-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT-NNNNN-GACTACHVGGGTATCTAATCC-3′). These primers include adapter sequences compatible with the Illumina library preparation workflow and are specific to the amplification of the bacterial V3–V4 regions. Polymerase chain reaction (PCR) was performed using template DNA (5 ng µl^−1^) with Blend Taq Plus polymerase (Toyobo, Osaka, Japan) following the manufacturer’s instructions. The PCR cycling conditions were initial denaturation at 94 °C for 2 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 30 s. The resulting PCR amplicons were verified by 1.5% agarose gel electrophoresis at 100 V for 25 min, stained with ethidium bromide, and visualized under UV illumination. The PCR amplicons were sent to Bioengineering Lab. Co., Ltd. (Kanagawa, Japan) for sequencing. The amplicons were purified using AMPure XP Beads (Beckman Coulter, Brea, CA, USA), and the DNA concentrations were measured using a Synergy H1 multimode microplate reader (Agilent Technologies, Santa Clara, CA, USA) and a QuantiFluor dsDNA System (Promega, Madison, WI, USA). Library preparation was performed using the purified amplicons with dual-index barcoding to enable sample multiplexing. After determining the library concentrations as described above, quality checking was performed using a dsDNA 915 Reagent kit (Agilent Technologies) on a Fragment Analyzer System (Agilent Technologies). Each library was pooled at an equimolar concentration and underwent 300 bp paired-end sequencing using MiSeq Reagent Kit v3 (Illumina, San Diego, CA, USA) on a MiSeq benchtop sequencer (Illumina).

Data analysis

To generate amplicon sequence variants (ASVs), DADA2 ver. 1.16 (Callahan et al., 2016) in RStudio software ver. 2023.12.0 + 369 was used for trimming and filtering the paired-end FASTQ reads data obtained following MiSeq high-throughput sequencing. Low-quality distributions and primers from each forward and reverse sequence were trimmed using parameters set to truncLen = c(290, 230) and trimLeft = c(17,21), respectively, and then the sequences were filtered using maxN = 0, maxEE = c(2,2), and truncQ = 2 parameters. Next, the forward and reverse reads were merged, after which chimeras and short reads (<400 bp) were discarded. The taxonomic classification of representative ASVs from phylum to genus level was assigned using the SILVA ver. 138.1 prokaryotic SSU database (Quast et al., 2013) as a reference dataset. Before downstream analysis, ASVs that were unclassified at the phylum level and were assigned as non-bacteria groups (e.g., chloroplasts and mitochondria) were manually removed.

To standardize sequencing depth at minimum read count across samples, abundance-based read resampling was performed using the rrarefy function of the vegan package ver. 2.6.4 (Oksanen et al., 2023) among all samples. Rarefaction curves were generated using the rarecurve function of the vegan package and the ggplot2 package ver. 3.4.2 (Wickham, 2016). The Coverage function in the entropart package ver. 1.6.12 (Marcon & Hérault, 2015) was used to calculate coverage and evaluate whether the sequencing depth was sufficient to fully represent the bacterial communities in each sample before and after rarefaction. After pooling ASVs at the lowest taxonomic level (bacterial genus), the bacterial community composition at the phylum and genus levels was visualized for each sample as bar plots using ggplot2. Taxa with a relative abundance of <1% across all samples were grouped and represented as “Others”.

Bacterial genera with a relative abundance >1% across all samples were defined as the core hindgut bacteria of Apis cerana japonica. All ASVs assigned to these core genera underwent similarity searches against the bacterial 16S rRNA gene sequence database in EzBioCloud (Yoon et al., 2017) to identify the closest related species.

Non-metric multidimensional scaling (NMDS) plot was generated based on Bray–Curtis dissimilarity using the metaMDS function from the vegan package and visualized using the ggplot2 package, in order to assess the β-diversity of the hindgut microbiota at the ASV level across the three time periods.

Statistical analysis

To detect differences in the hindgut microbial compositions among the three periods, we performed pairwise comparisons using permutational multivariate analysis of variance (PERMANOVA) based on the Bray–Curtis dissimilarity index with 9,999 permutations using the pairwise.adonis function of the pairwiseAdonis package ver. 0.4.1 (Martinez, 2020).

To investigate the effect of sampling period on the abundance of core bacterial genera, a generalized linear mixed model (GLMM) analysis was conducted, assuming a Poisson distribution with a log link function, using the glmer function from the lme4 package ver. 1.1.32 (Bates et al., 2015). The read count of each core bacterial genus was set as the response variable, with sampling period as a fixed effect and colony as a random effect. P-values < 0.05 were considered statistically significant for all comparisons.

Results

Sequence dataset overview

The high-throughput sequencing generated 1,381,026 raw reads (mean ± SD: 38,362 ± 5,199) from all samples. After trimming and filtering, 1,098,689 high-quality reads (mean ± SD: 30,519 ± 3,678) remained (Table S1), clustering into 260 ASVs (mean ± SD: 59 ± 10). Next, all nonbacterial reads were removed from all samples, and a read count-based cutoff was applied to match the minimum read count (21,808), resulting in a final dataset of 241 ASVs (mean ± SD: 54 ± 8). The rarefaction curves for each sample plateaued at the minimum read depth (Fig. S1), and the estimated coverage was exceeded 99% for all samples (Table S2), indicating that the sequencing depth was sufficient to identify most of the hindgut bacteria in the study samples.

Hindgut microbiota composition

The composition at the phylum and genus levels for each period were described for those taxa with a relative abundance >1% across all samples, while those with a relative abundance <1% were grouped as “Others” (Fig. 1). At the phylum level, the hindgut microbiota of A. c. japonica was dominated by Actinobacteriota (0.4–13.1%), Bacteroidota (11.7–36.3%), Firmicutes (9.2–33.9%), and Proteobacteria (35.2–68.6%), collectively accounting for >99.9% of the relative abundance. At the genus level, six bacterial genera, namely Apibacter (11.6–36.3%), Bifidobacterium (0.2–13.1%), Bombilactobacillus (0.3–13.6%), Gilliamella (26.6–59.8%), Lactobacillus (6.9–28.2%), and Snodgrassella (0.7–22.4%) predominated, accounting for >96% of the hindgut microbiota in all three periods and all samples. The relative abundance of “Unclassified”, of which sequences were not assigned at the genus level, and “Others” was 0.0–5.4% and 0.0–12.1%, respectively. The details of relative abundance at the phylum and genus levels are listed in Tables S3–S6.

Hindgut microbiota composition of the Japanese honey bees (Apis cerana japonica) sampled for each period at the (A) phylum and (B) genus levels.BO, before overwintering; OW, during overwintering; AO, after overwintering.

Further examination of the six major bacterial genera revealed that, except for genus Bifidobacterium, many ASVs assigned to the core genera exhibited sequence similarities with type strains in the EzBioCloud database (Yoon et al., 2017) below the 98.7% threshold, which is commonly used to distinguish closely related species (Chun et al., 2018) (Table 2). Notably, 90% (27/30) of the ASVs assigned to Gilliamella showed this pattern, followed by Snodgrassella (65%, 13/20), Bombilactobacillus (60%, 3/5), Apibacter (53.8%, 7/13), and Lactobacillus (41.6%, 5/12), suggesting potentially novel bacterial species within the hindgut microbiota of A. c. japonica.

Table 2: List of BLAST results against the EzBioCloud 16S rRNA database of ASVs assigned to the six core bacteria genera in the hindgut of the Japanese honey bees (Apis cerana japonica).

The result of the NMDS plot of β-diversity of the hindgut microbiota at the ASVs level based on Bray-Curtis dissimilarity is shown in Fig. 2. Pairwise comparisons among the three sampling periods revealed a significant difference in hindgut microbiota composition only between the BO and OW groups (F = 3.037, R^2^ = 0.121, p = 0.029; Table S7).

Nonmetric multidimensional scaling (NMDS) ordination plots of hindgut microbiota of the Japanese honey bees (Apis cerana japonica) at three sampling periods.The plot was generated with the Bray–Curtis dissimilarity index based on the ASVs obtained from each sample. BO, before overwintering; OW, during overwintering; and AO, after overwintering.

Comparison of the core genera among the three periods

The GLMM analysis revealed that the OW group had a significant positive effect on the read counts of Bifidobacterium, Bombilactobacillus, and Lactobacillus (coefficients: 0.977, 1.036, and 0.320; 95% CI [0.237–1.716], [0.138–1.933], and [0.131–0.509]; p = 0.009, 0.024, and 0.001, respectively; Table S8).

Data availability

The raw amplicon sequence datasets generated in this study are available in the DDBJ Sequence Read Archive (accession numbers: DRR685263–DRR685298 for DRA Run and PRJDB20791 for BioProject). All scripts and datasets are deposited to figshare under DOI: 10.6084/m9.figshare.29396408.

Discussion

This study revealed that the hindgut microbiota of Apis cerana japonica was dominated by four phyla: Actinobacteriota, Bacteroidota, Firmicutes, and Proteobacteria, and six core bacterial genera: Apibacter, Bifidobacterium, Bombilactobacillus, Gilliamella, Lactobacillus, and Snodgrassella. This finding is consistent with previous studies on the gut microbiota of honey bees (Kwong et al., 2017; Dong et al., 2020). In contrast, compared with the core gut microbiota of A. mellifera, Apibacter is more abundant in Asian honey bee species, such as A. cerana, A. dorsata, and A. andreniformis (Kwong & Moran, 2016; Kwong et al., 2017; Duong et al., 2020; Ellegaard et al., 2020; Khan et al., 2023). The hindgut microbiota of Apis cerana japonica showed a similar trend at the genus level, suggesting that host phylogeny influenced microbial community structure. However, the identification of ASVs with similarities lower than the threshold for distinguishing closely related species suggests the presence of many potentially novel bacterial species, despite their genus-level similarity. Further studies involving bacterial isolation, biochemical characterization, and genome analysis are warranted to elucidate the taxonomy and function of these candidate novel bacteria.

The hindgut microbiota composition of A. c. japonica in the OW group differed significantly from that of the BO group. Notably, the mean relative abundance of Bifidobacterium, Bombilactobacillus, and Lactobacillus in OW group was higher than that in BO group. These three core bacterial genera are known to produce SCFAs from pollen-derived polysaccharides and nectar-derived glucose (Zheng et al., 2017, 2019). Among the SCFAs derived from honey bee gut bacteria, butyrate is absorbed into the hemolymph via the ileum or rectum and is therefore considered an important energy source for thermogenesis to maintain hive temperature during the overwintering period (Den Besten et al., 2013; Zheng et al., 2017). Moreover, genera Bifidobacterium and Lactobacillus contribute substantially to infection control and immune regulation in honey bees through mechanisms such as antimicrobial activity against pathogens (Wu et al., 2013) and upregulating antimicrobial peptide expression (Daisley et al., 2020). Therefore, the genus-level increase in Bifidobacterium, Bombilactobacillus, and Lactobacillus in the OW group may play a beneficial role in the overwintering of honey bees in terms of thermogenesis and immune activation. This study has limitations, as it did not experimentally evaluate the SCFAs-producing capacities or the immunomodulatory effects of Bifidobacterium, Bombilactobacillus, and Lactobacillus. Further studies quantifying SCFA levels in the gut and analyzing the expression of immune-related genes during OW are necessary to clarify the functional roles of these gut bacteria in successful overwintering.

The observed compositional changes in the hindgut microbiota of Apis cerana japonica in the overwintering period are intriguing. A possible contributing factor is the difference in pollen and nectar sources consumed by honey bees before and during overwintering period. Honey bees forage across a wide temperature range (10–40 °C) (Abou-Shaara et al., 2017), but during the overwintering period, when temperatures fall below 10 °C, they rarely leave the hive to forage (Joshi & Joshi, 2010). Consequently, honey bees are more likely to consume stored pollen and honey during overwintering period. Furthermore, the consumption of aged or stored pollen and honey can influence gut microbiota composition through changes in physiological parameters (Maes et al., 2016). In our study, although daily maximum temperatures exceeded 10 °C on all sampling days before and after overwintering period, only one-third of the days during the overwintering period reached this threshold (Japan Meteorological Agency, 2025). Another factor that may influence hindgut microbiota is variations in hive temperature. Typically, the hive temperature is maintained at 33–35.5 °C (Abou-Shaara et al., 2017). In A. c. japonica, the average winter hive temperature is 30.7 °C, while the average temperature before and after winter is 34.3 °C (Akimoto, 2000). This temperature fluctuation may affect bacterial growth rates, thereby altering microbiota composition (Ludvigsen et al., 2015; Kešnerová et al., 2017).

It is important to note that slight differences in microbiota composition were observed among colonies. Differences in gut microbiota between colonies have been reported to be influenced by the diet collected from habitat-specific floral sources (Castelli et al., 2022; Ricigliano, Williams & Oliver, 2022) and host genotype (Bridson et al., 2022). Our findings highlight the necessity of considering variation among colonies when evaluating seasonal changes in the hindgut microbiota of honey bees.

Conclusions

This study on the hindgut microbiota of A. c. japonica revealed the influence of phylogeny on microbiota composition, the presence of potentially novel species, and distinctive compositional changes during the overwintering period. The biochemical properties of the genera that increased during overwintering period (i.e., genera Bifidobacterium, Bombilactobacillus, and Lactobacillus) suggest that these changes supply energy for thermogenesis and activate the host immune system. Further surveys in other regions with different dietary environments and studies focusing on elucidating the functional roles of hindgut microbiota during overwintering and their symbiotic relationship with host health are warranted.

Supplemental Information

10.7717/peerj.20050/supp-1Supplemental Information 1Supplementary figure and tables.

10.7717/peerj.20050/supp-2Supplemental Information 2Raw forward sequence data of X12-3.

10.7717/peerj.20050/supp-3Supplemental Information 3Raw forward sequence data of X3-1.

10.7717/peerj.20050/supp-4Supplemental Information 4Raw forward sequence data of X12-1.

10.7717/peerj.20050/supp-5Supplemental Information 5Raw reverse sequence data of X12-2.

10.7717/peerj.20050/supp-6Supplemental Information 6Raw forward sequence data of X10-1.

10.7717/peerj.20050/supp-7Supplemental Information 7Raw forward sequence data of X3-3.

10.7717/peerj.20050/supp-8Supplemental Information 8Raw reverse sequence data of X3-1.

10.7717/peerj.20050/supp-9Supplemental Information 9Raw reverse sequence data of X10-1.

10.7717/peerj.20050/supp-10Supplemental Information 10Raw reverse sequence data of X12-1.

10.7717/peerj.20050/supp-11Supplemental Information 11Raw forward sequence data of X3-2.

10.7717/peerj.20050/supp-12Supplemental Information 12Raw forward sequence data of X10-3.

10.7717/peerj.20050/supp-13Supplemental Information 13Raw reverse sequence data of X10-3.

10.7717/peerj.20050/supp-14Supplemental Information 14Raw reverse sequence data of X10-2.

10.7717/peerj.20050/supp-15Supplemental Information 15Raw reverse sequence data of X3-2.

10.7717/peerj.20050/supp-16Supplemental Information 16Raw reverse sequence data of X3-3.

10.7717/peerj.20050/supp-17Supplemental Information 17Raw forward sequence data of X12-2.

10.7717/peerj.20050/supp-18Supplemental Information 18Raw forward sequence data of X10-2.

10.7717/peerj.20050/supp-19Supplemental Information 19Raw reverse sequence data of X12-3.

10.7717/peerj.20050/supp-20Supplemental Information 20Raw reverse sequence data of T10-2.

10.7717/peerj.20050/supp-21Supplemental Information 21Raw reverse sequence data of T3-2.

10.7717/peerj.20050/supp-22Supplemental Information 22Raw reverse sequence data of T12-1.

10.7717/peerj.20050/supp-23Supplemental Information 23Raw reverse sequence data of T3-1.

10.7717/peerj.20050/supp-24Supplemental Information 24Raw reverse sequence data of T3-3.

10.7717/peerj.20050/supp-25Supplemental Information 25Raw reverse sequence data of T10-1.

10.7717/peerj.20050/supp-26Supplemental Information 26Raw forward sequence data of T10-1.

10.7717/peerj.20050/supp-27Supplemental Information 27Raw forward sequence data of T10-2.

10.7717/peerj.20050/supp-28Supplemental Information 28Raw forward sequence data of T12-1.

10.7717/peerj.20050/supp-29Supplemental Information 29Raw forward sequence data of T12-2.

10.7717/peerj.20050/supp-30Supplemental Information 30Raw forward sequence data of T3-2.

10.7717/peerj.20050/supp-31Supplemental Information 31Raw forward sequence data of T10-3.

10.7717/peerj.20050/supp-32Supplemental Information 32Raw reverse sequence data of T12-2.

10.7717/peerj.20050/supp-33Supplemental Information 33Raw reverse sequence data of T12-3.

10.7717/peerj.20050/supp-34Supplemental Information 34Raw reverse sequence data of T10-3.

10.7717/peerj.20050/supp-35Supplemental Information 35Raw forward sequence data of T12-3.

10.7717/peerj.20050/supp-36Supplemental Information 36Raw forward sequence data of T3-1.

10.7717/peerj.20050/supp-37Supplemental Information 37Raw forward sequence data of T3-3.

10.7717/peerj.20050/supp-38Supplemental Information 38Raw forward sequence data of I3-1.

10.7717/peerj.20050/supp-39Supplemental Information 39Raw reverse sequence data of I12-3.

10.7717/peerj.20050/supp-40Supplemental Information 40Raw forward sequence data of I12-3.

10.7717/peerj.20050/supp-41Supplemental Information 41Raw forward sequence data of I3-2.

10.7717/peerj.20050/supp-42Supplemental Information 42Raw reverse sequence data of I12-1.

10.7717/peerj.20050/supp-43Supplemental Information 43Raw reverse sequence data of I10-3.

10.7717/peerj.20050/supp-44Supplemental Information 44Raw reverse sequence data of I10-1.

10.7717/peerj.20050/supp-45Supplemental Information 45Raw forward sequence data of I3-3.

10.7717/peerj.20050/supp-46Supplemental Information 46Raw reverse sequence data of I3-2.

10.7717/peerj.20050/supp-47Supplemental Information 47Raw forward sequence data of I10-3.

10.7717/peerj.20050/supp-48Supplemental Information 48Raw reverse sequence data of I3-1.

10.7717/peerj.20050/supp-49Supplemental Information 49Raw forward sequence data of I10-1.

10.7717/peerj.20050/supp-50Supplemental Information 50Raw reverse sequence data of I3-3.

10.7717/peerj.20050/supp-51Supplemental Information 51Raw forward sequence data of I12-2.

10.7717/peerj.20050/supp-52Supplemental Information 52Raw reverse sequence data of I12-2.

10.7717/peerj.20050/supp-53Supplemental Information 53Raw forward sequence data of I10-2.

10.7717/peerj.20050/supp-54Supplemental Information 54Raw reverse sequence data of I10-2.

10.7717/peerj.20050/supp-55Supplemental Information 55Raw forward sequence data of I12-1.

10.7717/peerj.20050/supp-56Supplemental Information 56Raw forward sequence data of H10-1.

10.7717/peerj.20050/supp-57Supplemental Information 57Raw forward sequence data of H10-2.

10.7717/peerj.20050/supp-58Supplemental Information 58Raw forward sequence data of H12-2.

10.7717/peerj.20050/supp-59Supplemental Information 59Raw reverse sequence data of H3-2.

10.7717/peerj.20050/supp-60Supplemental Information 60Raw reverse sequence data of H12-1.

10.7717/peerj.20050/supp-61Supplemental Information 61Raw reverse sequence data of H12-3.

10.7717/peerj.20050/supp-62Supplemental Information 62Raw reverse sequence data of H10-2.

10.7717/peerj.20050/supp-63Supplemental Information 63Raw forward sequence data of H3-3.

10.7717/peerj.20050/supp-64Supplemental Information 64Raw forward sequence data of H10-3.

10.7717/peerj.20050/supp-65Supplemental Information 65Raw forward sequence data of H3-1.

10.7717/peerj.20050/supp-66Supplemental Information 66Raw forward sequence data of H12-1.

10.7717/peerj.20050/supp-67Supplemental Information 67Raw reverse sequence data of H10-1.

10.7717/peerj.20050/supp-68Supplemental Information 68Raw reverse sequence data of H10-3.

10.7717/peerj.20050/supp-69Supplemental Information 69Raw reverse sequence data of H3-3.

10.7717/peerj.20050/supp-70Supplemental Information 70Raw reverse sequence data of H12-2.

10.7717/peerj.20050/supp-71Supplemental Information 71Raw forward sequence data of H3-2.

10.7717/peerj.20050/supp-72Supplemental Information 72Raw forward sequence data of H12-3.

10.7717/peerj.20050/supp-73Supplemental Information 73Raw reverse sequence data of H3-1.