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

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.

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⁻¹ Tris-HCl and 1 mmol L⁻¹ 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⁻¹) 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.

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.

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.

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² = 0.121, p = 0.029; Table S7).

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.