The adaptation of bumblebees to extremely high elevation associated with their gut microbiota

ABSTRACT Bumblebees are among the most abundant and important pollinators for sub-alpine and alpine flowering plant species in the Northern Hemisphere, but little is known about their adaptations to high elevations. In this article, we focused on two bumblebee species, Bombus friseanus and Bombus prshewalskyi, and their respective gut microbiota. The two species, distributed through the Hengduan Mountains of southwestern China, show species replacement at different elevations. We performed genome sequencing based on 20 worker bee samples of each species. Applying evolutionary population genetics and metagenomic approaches, we detected genes under selection and analyzed functional pathways between bumblebees and their gut microbes. We found clear genetic differentiation between the two host species and significant differences in their microbiota. Species replacement occurred in both hosts and their bacteria (Snodgrassella) with an increase in elevation. These extremely high-elevation bumblebees show evidence of positive selection related to diverse biological processes. Positively selected genes involved in host immune systems probably contributed to gut microbiota changes, while the butyrate generated by gut microbiota may influence both host energy metabolism and immune systems. This suggests a close association between the genomes of the host species and their microbiomes based on some degree of natural selection. IMPORTANCE Two closely related and dominant bumblebee species, distributed at different elevations through the Hengduan Mountains of southwestern China, showed a clear genomic signature of adaptation to elevation at the molecular level and significant differences in their respective microbiota. Species replacement occurred in both hosts and their bacteria (Snodgrassella) with an increase in elevation. Bumblebees’ adaptations to higher elevations are closely associated with their gut microbiota through two biological processes: energy metabolism and immune response. Information allowing us to understand the adaptive mechanisms of species to extreme conditions is implicit if we are to conserve them as their environments change.


Metagenomic sequencing of bumblebees
Deep-coverage metagenomic shotgun sequencing performed on all 20 samples for each species (Fig. 1; Table S1) resulted in an average of 36 Gbp for each sample.Subsequent processing divided each sample datum into two parts (Table S2).One part derived from the bumblebee host showed an average of 28.36 Gbp, and the second part derived from gut microbes had an average of 8.39 Gbp.Among 40 samples, 35 (N = 19 for B. friseanus; N = 16 for B. prshewalskyi) had over 15 Gbp corresponding to the host.The remaining five had 3. 45-11.96Gbp corresponding to the host.All 40 samples had over 1.2 Gbp corresponding to the gut microbiome.

Host population structure and adaptation
Genome alignment indicated an average of 112.5× sequencing coverage for each individual relative to a 225-Mb reference genome of Bombus pyrosoma.The average nucleotide identity (ANI) between B. friseanus and B. prshewalskyi genomes was 99.03%.A total of 112,132 single nucleotide polymorphisms (SNPs) were identified.Reduced genome-wide linkage disequilibrium and genetic diversity (π) were observed in the B. prshewalskyi populations at extremely high elevations (Fig. S1A and B).Based on linkage disequilibrium, the estimated effective population sizes of B. friseanus and B. prshewalskyi populations were 1,065 and 69, respectively.Population genetic analyses showed that B. friseanus and B. prshewalskyi populations formed two branches and divided distinctly in the maximum likelihood tree.This clear division was demonstrated further by genetic structure analysis at K = 2 (Fig. 2A).Principal component analysis (PCA) also supported the clear division of these two species.In PCA, the first principal component (PC) axes explained 11.18% of the total variation separating both bumblebee populations (Fig. S1C).

Microbiota composition and functions
Taxonomic profiles of bumblebee gut microbiomes estimated using the Kraken2 software showed that the core gut bacteria of bumblebees (i.e., Gilliamella, Snodgrassella, Lactobacillus, Bifidobacterium, and Apibacter) dominated most samples, although seven samples were dominated exclusively by pathogens.Specifically, four samples from B. friseanus were dominated by eukaryotic Crithidia sp., and three samples of B. prshewalskyi were dominated by potentially pathogenic bacteria including Pseudomonas sp., Serratia sp., and Rahnella sp. (Fig. 3A).
PCoA based on the distance calculated by KmerFreqCalc showed distinct clusters of the gut microbiomes from two bumblebee species.The first PC explained 43.5% and the second 34.1% (Fig. 3B).The two clusters were significantly associated with host specificity [permutational analysis of variance (PERMANOVA) F 1, 29 = 14.30, r 2 = 0.33, P < 0.01).In addition, the difference between gut microbiomes showed a weak relation to the plant species foraged on by the worker castes of both species (PERMANOVA F 8, 22 = 3.07, r 2 = 0.53, P < 0.05).

Metagenomic binning
Metagenomic binning yielded 23 bins showing >85% completeness and <5% contami nation (Table S5).Taxonomic identification classified eight bins representing genera of eusocial bee gut core bacteria.This represented one bin from Gilliamella, two from Snodgrassella, two from Lactobacillus, two from Bifidobacterium, and one from Apibacter.Phylogenetic analysis showed that most of the core bacterial bins clustered together within typical strains of bumblebee gut bacteria, forming distinct clusters with their relatives in the honeybee gut, excluding the bin of Snodgrassella_1 located at the base of the Snodgrassella clade (Fig. S2).Quantitative analysis showed that the relative abun dance patterns of four bins (Snodgrassella_1, Snodgrassella_2, Bifidobacterium_2, and Lactobacillus_2) differed significantly in the two bumblebee species (Wilcoxon rank-sum test, P < 0.01) (Fig. 4A; Fig. S3).Notably, Snodgrassella bins showed clear host specificity.Snodgrassella_1 dominated B. friseanus, while Snodgrassella_2 dominated B. prshewalskyi.

Two Snodgrassella bins
The completeness of Snodgrassella_1 and Snodgrassella_2 bins was 94.77% and 88.82%, respectively.The ANI between the two bins was 77.65%.Phylogenomic analysis indica ted that the two Snodgrassella bins differed significantly.Snodgrassella_1 clustered on the base of the Snodgrassella clade, while Snodgrassella_2 clustered with other strains from bumblebees (Fig. S2).Alignment with the genome of Snodgrassella alvi_wkB2 showed that each bin had unique genes (59 for Snodgrassella_1 and 336 for Snodgras sella_2), with 1,377 genes shared (Fig. 4B; Table S6).Five pathways showed enrichment in the 336 Snodgrassella_2 unique genes, including metabolic pathways, microbial metabolism in a diverse environment, butanoate metabolism, oxidative phosphoryla tion, and sulfur metabolism.No pathways were enriched in the 59 genes unique to Snodgrassella_1 (Fig. 4C).Further analysis indicated that genes involved in butanoate metabolism showed significant differences in the two bins (Fig. S4A).Snodgrassella_2 had twice as many genes involved in the butanoate synthesis pathway as Snodgras sella_1 (Fig. S4B).
To identify the adaptation signals of the two Snodgrassella bins, 15 sets of sequencing reads for each of the Snodgrassella bins were derived from raw data, identifying a total of 95,215 SNPs.Population genetic analyses, including population genetic structure (K = 2) and the maximum likelihood tree, showed a clear division between Snodgrassella_1 and Snodgrassella_2 (Fig. 4D).The MK test identified 22 candidate genes under positive selection from 496 with sufficient polymorphisms (Table S7; Fig. 4E).These genes involved in multiple KEGG pathways, including beta-lactam resistance (ko01501: two genes), two-component systems (ko02020: four genes), ABC transporters (ko02010: two genes), purine metabolism (ko00230: two genes), amino acid metabolism (ko00250: one gene, ko00260: one gene, and ko00400: one gene), riboflavin metabolism (ko00740: one gene), pentose phosphate pathway (ko00030: one gene), and homologous recombina tion (ko03440: one gene).

DISCUSSION
We focused on two bumblebee species from the same subgenus as well as their gut microbiota, allowing us to hypothesize that genetic diversity in these two prospective host insects occurs along distinct zoogeographic ranges as shaped by adaptive selection.We also expected to detect differences between their gut microbiotas via metagenomic approaches, which may be involved in host genome differences.
We identified 23 genes under positive selection in B. prshewalskyi living at extremely high elevations.Consequently, it is not surprising that the gut microbiotas of both bumblebee species are significantly different in composition and function.Under positive selection, these 23 genes in B. prshewalskyi underlie diverse biological pro cesses, implying that bumblebees adapt to extremely high elevations via complex and systematic mechanisms.Three genes associated with the regulation of cellular respira tion (TRAP1 and PNPT1) and response to stress/oxidative stress (TRAP1, PNPT1, and PERO) significantly reflect the adaptation to hypoxia at greater heights.In particular, the positively selected BRM gene indicates that transcriptional regulation plays an important role in maintaining normal development and growth of bumblebees living under harsh conditions (i.e., lower ambient temperatures and hypoxia).A recent study on males of B. terrestris in Britain also showed that environmental pressures most likely contributed to recent changes in genes underlying physiology, neurology, and wing development (9).
Furthermore, BRM together with Spätzle (Spz) suggests ongoing stress on the bumblebees' immune systems, with BRM regulating the innate immune response and SPZ in the Toll-like signaling pathway.The Toll-like signaling pathway is one of the most important intracellular signaling pathways in innate immunity systems in bees (24,25) and other insects including Drosophila (26).In humans, acute high-elevation exposure upregulated the inflammatory signaling pathways and may have sensitized the toll-like receptor 4 (TLR4) signaling pathway to subsequent inflammatory stimuli (27).Further evaluation of the expression level and function of BRM or Spz and other immunological tests may provide additional support via future experimentation.Seven genes (FAS, FASC, OGFD1, LAC1, INSR, GRB14, and Zn236) involved in metabolic processes imply strong pressures on energy metabolism in montane bumblebees.Among them, INSR (insulin-like peptide receptor) and GRB14 are mentioned because they are involved in the insulin or insulin-like growth peptide signaling (IIS) pathway, an evolutionar ily conservative nutrient-sensing pathway that modulates energy metabolism and development in metazoans (28,29).In honeybees, the IIS pathway participates in the regulation of caste development (30)(31)(32) and the response to cold stress (33).The IIS pathway is also reported to contribute to the adaptive response to hypoxia in other insects including Drosophila (34,35) and Tibetan locusts (6).Consequently, we propose here that the genetic variants of INSR and GRB14 in B. prshewalskyi play important roles in adapting bumblebee workers to hypoxia at extremely high elevations by regulating energy metabolism.
Additionally, the results suggest a potential correlation between the gut microbiota and those bumblebee adaptations, focusing on butanoate metabolism in the microbiota of B. prshewalskyi and the bins of Snodgrassella_2.Butanoate (or butyrate) is one of the short-chain fatty acids (SCFAs) generated by gut microbes.It can activate G-cou pled-receptors directly, inhibit histone deacetylases, and serve as an energy substrate, regulating an animal's physiology and health (36)(37)(38).In humans, butyrate is the primary energy source for colonocytes and also protects against colorectal cancer and inflam mation (39,40).It can improve the host's metabolism of carbohydrates and lipids by serving as a signal molecule (41).Mice fed with a butyrate-enriched high-fat diet showed increased thermogenesis and energy expenditure and appeared resistant to obesity (42).Though simpler than the human or mice gut microbiota, the honeybee gut microbiota exhibits a variety of beneficial effects, from gut-centric functions like digestion (43), detoxification (44), and defense from pathogens (45), to more peripheral processes like behavior (46).The bumblebee gut microbiota has also been shown to provide some resistance to a trypanosomatid parasite (Crithidia bombi) (47).These effects are thought to be mediated, in part, by SCFAs (43).Considered synthetically, we propose that butyrate concentrations in the guts of bumblebees affect their energy metabolism and gut immunity.As both processes appear to be under stress at extremely high elevations, the results suggest that the gut microbiota may influence their hosts' adaptations by adjusting the concentration of butyrate.Biochemical details of the impact effect of the gut microbiota on the bumblebee host should be elucidated further by experimentation on the gut butyrate and host physiology of these two bumblebee species or on related model species.
The different gut communities in the two bumblebee species reflect the hosts' distinctive influence on their own gut microbiota.It is well known that host genetic background and diet are major factors influencing the gut microbiota (48)(49)(50)(51).As plant species foraged upon by the bumblebee, samples were shown to correlate weakly with differences in gut microbiotas, and these differences may be mainly due to the hosts' genetic background changes.The 22 genes involved in multiple KEGG pathways under positive selection in Snodgrassella_2 suggest strong host stress on their signaling pathway and antibiotic resistance process.Snodgrassella bacteria, forming a biofilm on the bumblebees' ileum wall (15), directly encounter their host's immune system.A recent study reported that the host specificity of honeybee gut bacteria was determined through reactive oxygen species that are regulated by immune deficiency and Toll pathways (52).Therefore, those host genetic variants involved in the immune system probably contribute to the species replacement of Snodgrassella.
In conclusion, we have helped to show that at extremely high elevations, B. prshewal skyi shows signs of positive selection in diverse biological processes, including cellular respiration, stress response, growth and development, locomotory behavior, immune system, and energy metabolism.We also found differences in both the composition and function of the bumblebee gut microbiome along an elevation gradient.Differen ces between the host bee species and their gut microbe species, as driven by adap tive selection, are probably functional and suggest a closer association between the genomes of hosts and their microbiomes.As adaptations shared by these bumblebee species with their microbiota indicate some degree of specialization, it may also reflect some degree of coevolution, but this will require further study.

Bumblebee sampling
We  S1 for the detailed information of sampling sites).As bumblebees are eusocial insects, we tried to avoid collecting individuals from the same colony by sampling more than one site on the three mountains for both species.The subsequent population genetic analyses of both species showed a high divergence among individual workers, suggesting that we did collect specimens from more than one colony (Fig. 2A).
We recorded the observation date, time, and plant species on which each foraging worker was collected.Each specimen was netted and stored in a separate 50 mL Eppendorf tube and euthanized at 4°C prior to gut dissection.Specimens were identified as species using morphological characters (18) under field conditions.Gut dissections and their preservation were performed in the field on the day of collection using a Stereo Microscope (Zeiss Stemi 508) as follows.
The surface of each bumblebee was first washed in droplets of 70% ethanol for 30 s and then rinsed with sterile water three times to remove external contaminants.We then pulled out the entire gut from the abdomen terminus using sterilized forceps.Each gut specimen was rinsed twice immediately with 0.9% sterile NaCl solution to maintain cell pressure.Each sample was stored separately in a 2mL cryotube placed in liquid nitrogen.After returning to the laboratory, gut specimens were frozen at -80°C until DNA extraction.In total, we collected 20 gut samples for each Bombus species.

Metagenomic shotgun sequencing and processing
Total genomic DNA was extracted from gut materials using the E.Z.N.A. Soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA).The concentration and purity of extracted DNA were determined with TBS-380 and NanoDrop2000, respectively.DNA extract quality was checked on a 1% agarose gel.The DNA extract was fragmented to an average size of about 400 bp using Covaris M220 (Gene Company Limited, China) for paired-end library construction.The paired-end library was constructed using NEXTFLEX Rapid DNA-Seq (Bioo Scientific, Austin, TX, USA).Adapters containing the full complement of sequenc ing primer hybridization sites were ligated to the blunt end of fragments.Paired-end sequencing was performed on Illumina NovaSeq 6000 (Illumina Inc., San Diego, CA, USA) at Majorbio Bio-Pharm Technology Co., Ltd.(Shanghai, China) using NovaSeq Reagent Kits, generating an average of 36 gigabyte bases.
Reads were first trimmed off with Trimmomatic v0.39 (parameters used: LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:40) (53) and inspected with FastQC v0.11.8 (54) to validate the trimming quality.Clean reads of each metagenomic sample were mapped against the reference genome of B. pyrosoma (GCA_14825855.1) using BWA v0.7.17-r1188 (55) to separate reads from the host and its microbes.B. pyrosoma belongs to the same subgenus as the two species, and its genome was annotated using Trinotate (http://trinotate.sourceforge.net/).About 80% of the total reads were derived from the host, and most of the remaining reads belonged to gut microbes.These two datasets were used in the subsequent bumblebee genomic and gut microbiome analyses.
PCA for the first two components was performed based on the SNP site of the bumblebee in Plink (62).The resulting graph was plotted using the ggplot2 package in R (63).
The population structure was determined using ADMIXTURE v 1.3.0(64).The estimation was K = 2 to K = 4, and the cross-validation (CV) error estimation was minimized at K = 2.A maximum-likelihood phylogenetic tree was constructed by the RAxML software (65), with the GTRGAMMA model and 100 bootstraps.An ascertainment bias correction was performed to correct for the impact of invariable sites in the data.
The MK test (8) was used to make inferences regarding historical selection by comparing fixed divergence and polymorphism at synonymous and nonsynonymous sites using a 2 × 2 contingency table.Significance was tested using Fisher's exact test.To correct the false positives, the proportion of amino acid fixations driven by positive selection (α) was calculated following the procedure proposed by Shapiro et al. (66).

FIG 2
FIG 2 Host population differentiation and adaptation to higher (4,200 m) elevations.(A) Maximum likelihood tree and population structure of two bumblebee species.(B) The result of the McDonald-Kreitman test on coding genes of bumble bees.Orange diamonds indicate those genes under positive selection, while gray circles indicate other genes.F N / F S : the fixed nonsynonymous/synonymous ratio between species; P N / P S : the polymorphism nonsynonymous/synonymous ratio within species.(C) Candidate genes under positive selection are involved in diverse gene ontology (GO) processes.Lines indicate that genes are associated with the GO processes.

FIG 3
FIG 3 Gut communities of two bumblebee species differing in composition and function.(A) Taxonomic profiles of bumblebee gut microbiota estimated using the Kraken2 software.(B) Principal coordinate analysis (PCoA) based on the distance calculated by KmerFreqCalc.The clusters of bumblebee gut microbiotas were significantly associated with host specificity [permutational analysis of variance (PERMANOVA) F 1, 29 = 14.30, r 2 = 0.33, P < 0.01] and related weakly with plant species foraged on by worker bees (PERMANOVA F 8, 22 = 3.07, r 2 = 0.53, P < 0.05).(C) Kyoto Encyclopedia of Genes and Genomes pathways with differential relative abundances between gut microbiomes of B. friseanus and B. prshewalskyi.

FIG 4
FIG 4 Metagenome binning revealed a pair of vicariously distributed species in Snodgrassella bacteria and the adaptation of Snodgrassella bacteria to the gut environment in a host restricted to a higher elevation.(A) The abundance of metagenome assembly bins in each healthy sample of bumblebees.The bins with significantly different abundances in two bumblebee species are marked with ** (P < 0.01, rank-sum test) or *** (P < 0.001, rank-sum test).(B) Comparison of gene compositions between two Snodgrassella bins using the genome of Snodgrassella alvi_wkB2 as the reference.(C) KEGG pathway enrichment analysis in the unique genes of Snodgrassella bins.(D) Maximum likelihood tree and population structure of two Snodg rassella bins.(E) The result of the McDonald-Kreitman test on genes in Snodgrassella bins.F N / F S : fixed nonsynonymous/synonymous ratio between species; P N / P S : polymorphism nonsynonymous/synonymous ratio within species.Colored diamonds indicate the genes under positive selection ((F N / F S ) / (P N / P S ) >1 and Fisher's exact test P < 0.05) in diverse KEGG pathways, while gray circles indicate other genes.
collected workers of B. friseanus and B. prshewalskyi for gut extraction from 11 August to 15 August 2020.Workers of B. friseanus were collected between 2,700 and 3,200 m on Yulong Snow Mountain and the Shika Snow Mountain, northwest Yunnan.Collections of B. prshewalskyi came from 4,200 to 4,700 m on Baima Snow Mountain (Fig. 1, and Table