Genetic Modification of a Hox Locus Drives Mimetic Color Pattern Variation in a Highly Polymorphic Bumble Bee

Abstract

Müllerian mimicry provides natural replicates ideal for exploring mechanisms underlying adaptive phenotypic divergence and convergence, yet the genetic mechanisms underlying mimetic variation remain largely unknown. The current study investigates the genetic basis of mimetic color pattern variation in a highly polymorphic bumble bee, Bombus breviceps (Hymenoptera, Apidae). In South Asia, this species and multiple comimetic species converge onto local Müllerian mimicry patterns by shifting the abdominal setal color from orange to black. Genetic crossing between the orange and black phenotypes suggested the color dimorphism being controlled by a single Mendelian locus, with the orange allele being dominant over black. Genome-wide association suggests that a locus at the intergenic region between 2 abdominal fate-determining Hox genes, abd-A and Abd-B, is associated with the color change. This locus is therefore in the same intergenic region but not the same exact locus as found to drive red black midabdominal variation in a distantly related bumble bee species, Bombus melanopygus. Gene expression analysis and RNA interferences suggest that differential expression of an intergenic long noncoding RNA between abd-A and Abd-B at the onset setal color differentiation may drive the orange black color variation by causing a homeotic shift late in development. Analysis of this same color locus in comimetic species reveals no sequence association with the same color shift, suggesting that mimetic convergence is achieved through distinct genetic routes. Our study establishes Hox regions as genomic hotspots for color pattern evolution in bumble bees and demonstrates how pleiotropic developmental loci can drive adaptive radiations in nature.

Keywords: Hox; Müllerian mimicry; bumble bees; genomic hotspot; noncoding RNA.

Fig. 1.

Mimetic color pattern diversity of the South Asia bumble bee mimicry group. A) Phylogenic position and color pattern diversity of the comimetic and distantly related B. breviceps (B. B), B. haemorrhoidalis (B. H), and B. trifasciatus (B. T) lineages, belonging to subgenera Alpigenobombus, Orientalibombus, and Megabombus, respectively. Note that this mimetic group displays parallel color variation on different body segments including thoracic, pleuron, anterior metasomal tergal (T1-3), and finally, 5th metasomal tergal (T5, also known as “tail” in bumble bees) regions, which displays the orange/black dimorphism studied here. All mimetic color patterns of these lineages are divided into 2 groups based on color phenotype (orange or black) on metasomal T5, highlighted by orange and gray shading, respectively. The color pattern diagram with a red star represents B. grahami (B. G.), the sister species of B. breviceps. B) Major mimetic color patterns of the 3 lineages and their simplified geographic distribution. Dots and connecting numbers on the map indicate different locations where mimicry among the 3 lineages cooccur. Graphs below the maps illustrate major mimetic color patterns in those respective locations by comimetic species.

Fig. 2.

The Southern China orange/black color pattern transition zone of B. breviceps. Dots and squares on the map indicate localities from where bees were collected and phenotyped. The color pattern transition zone is highlighted by dashed white circles with distances labeled. Note that this transition zone is ∼300 km in width from northwest to southeast and 600 km in length from northeast to southwest, involving southern Yunnan, all of Guizhou, and southern Hunan. This transition zone occurs across altitudinal gradients, with the orange/black ratio gradually decreasing from northwest to southeast, largely concordant with the topographic transition from the Yun-Gui highlands to the lower elevation mountainous area of southeastern China. Orange and black color indicates orange- and black-tailed patterns, respectively. Pie charts indicate frequency of orange- and black-tailed patterns at each location, represented by percentage of orange-/black-tailed bees in the total collection from that location (supplementary data, Supplementary Material online). The size and portion of these charts reflect the number of samples per locality. Dots with a yellow circle indicate locations where samples were collected for GWAS and genetic crossing. The orange black transitions occur across an altitudinal gradient. Across the collection regions, each of the orange- and black-tailed patterns display variations on yellow banding numbers on anterior abdominal terga (T1-2), based on which each of them can be further partitioned into 1 yellow-banded and 2 yellow-banded color pattern groups, indicated by dots and squares, respectively. However, only the 1 yellow-banded pattern (indicated by red stars) of the orange/black phenotypes occur in the transition zone and was sampled for GWAS. Dots with yellow circles indicate sites where queens of the offspring males used for GWAS were collected. MM, Myanmar; THA, Thailand; LAO, Laos; VN, Vietnam; TB, Tibet Plateau; SC, Sichuan Basin; YG, Yunnan-Guizhou Plateau; NL, Nanling Mountains; WY, Wuyi Mountains.

Fig. 3.

The color locus revealed by GWAS. The Manhattan plot on the bottom is the whole-genome association showing that sites of fixed association (indicated in red) are confined to a single locus on chromosome 18. The plot above shows association within this chromosome, showing that fixation is confined to a ∼17.4 kb block 3′ of Abd-B (highlighted gray). The region highlighted in red indicates the position of the 1.4 kb nonrecombinant region that remains fixed with additional genotyping. The region highlighted in blue indicates the color locus controlling the red/black midabdominal color dimorphism in B. melanopygus (Tian et al. 2019). SNPs (black bars) and indels (orange bars) and their position are indicated below the 1.4 kb locus. Note that the black pattern contains 3 and 46 bp deletions at the nonrecombinant locus. The sequences show genotypes of different orange and black color patterns at the 1.4 kb locus, as well as that of B. grahami (B.G.) with only partial fragments containing the fixed SNPs and indels shown. Stars indicate patterns used for GWAS. Dots indicate color patterns used for narrowing the locus. Triangles indicate additional geographic color pattern genotypes for validating the association of this 1.4 kb locus with color morphs. Note that the locus remains fixed across all orange- and black-tailed patterns regardless of their geographic locations and variations on other segmental phenotypes (e.g. thorax and pleuron).

Fig. 4.

Phylogeographic relationships of different color pattern lineages of B. breviceps. A) Phylogenetic clustering of color pattern at the 1.4 kb color locus. Trees were constructed using ML methods with terminal nodes colored by the 2 tail color phenotypes. Triangular clades represent shared haplotypes, with size reflecting the frequency of the haplotype. Colored bars indicate broad geographic regions. B) Haplotype networks of the color locus. Each node represents a haplotype, lines between nodes represent a single base change, and the size of the node represents haplotype frequency. Nodes are colored by the color patterns containing that haplotype. Geographic regions are indicated using dashed polygons.

Fig. 5.

Gene expression analysis for the bithorax complex (BX-C) during bumble bee development. A) Temporal expression pattern of the abd-A and Abd-B genes at the dimorphic metasomal tergal segment 5 (T5) across pupal and early adult development. Images for pupal (P4-P15) and adult (QA-C12) stages used for qPCR analysis are provided below the gene expression plot. These images also demonstrate the process of setal pigmentation. The white arrows indicate polymorphic metasomal T5. Note that the orange and black setae are distinguishable starting from QA stage. At the QA stage, black setae exhibit a grayer tinge, while the orange setae become more golden orange. The 2 phenotypes become more distinguishable at the time of eclosion (0 h callow). Color continues to intensify in the first 24 h posteclosion. Both orange and black setae attain most of their color within the first 6 h of posteclosion and complete pigmentation by 12–24 h, largely in line with what has been found in other bumble bees (Tian and Hines 2018). B) The gene transcription profile of orange- and black-tailed bees across the BX-C complex at the QA period. The numbered line in the middle represents position of the chromosomal 18 corresponding to the BX-C complex. The red bar indicates the color locus. Positions of the infra-abdominal (iab) domain are labeled above the line. Gray dashed lines immediately below the chromosomal line indicate position of the iab8 lncRNA and mir-iab4 microRNA discovered by previous studies. Black solid lines indicate position of the protein coding sequence of BX-C Hox genes, as well as ncRNAs recognized by current studies. The iab6 lncRNA were shown in red solid line. The bottom represents transcriptome sequencing reads that mapped to the BX-C regions, with the height of the peak indicating read coverage. Each line represents the transcription profile of an individual bee, with orange- and black-colored graphs representing orange-tailed and black-tailed bees, respectively. Read coverage of the iab6 ncRNA is highlighted by a dashed square. Note that this lncRNA only has read coverage in the orange-tailed bees but little to no read coverage in the black-tailed samples, suggesting its tight association with the orange pigmentation. C) Gene expression of iab6 lncRNA in the dimorphic T5 and monomorphic metasomal T3. D) Heatmap indicating iab6 lncRNA expression level across abdominal tergites and pupal stages of B. breviceps. Note that this lncRNA is highly stage- and tissue-specific, with peak expression at T5 at the QA stage. Error bars in (A) and (C) indicate standard error (SE) of the mean. Stars above the error bars indicate levels of statistical significance *P < 0.05, **P < 0.01, ***P < 0.001 by 2-sample t-test.

Fig. 6.

Genetics and inheritance of orange/black dimorphism in B. breviceps comimics. A) Allelic variants in the color locus of B. breviceps (BB) and the comimetic B. montivagus (BM) and B. haemorrhoidalis (BH). Orange and gray shading indicates SNPs of the orange- and black-tailed form, respectively. Blue shading indicates alternative SNPs to those in B. breviceps. In B. montivagus, no fixed SNP or indels are similarly differentiated by color forms. In B. haemorrhoidalis, only black-tailed phenotypes were sampled and genotyped, and only a few sites of these black formed bees share the same alleles with the black-tailed B. breviceps. B) The left graph shows color patterns of males produced by 2 black-tailed B. montivagus workers. One black-tailed B. montivagus worker produced orange-tailed males, which is not consistent with the B. breviceps inheritance model in which black tails are controlled by recessive alleles and black workers are homozygous for rr and can produce only black r males. This indicates that the orange/black dimorphism in metasomal T5 is governed by different inheritance rules in B. montivagus.

Fig. 7.

The proposed mode of genetic control of the orange/black color dimorphism in B. breviceps. The schematic diagram illustrates how mutations in the color locus may drive orange/black setal pigmentation changes. This model proposes that the color locus of B. breviceps may correspond to an iab genetic switch that controls iab6 lncRNA expression in adult epidermis of metasomal T5. A mutation on the locus may cause the iab6 lncRNA switch to turn off in the black-tailed form, which causes expression silencing of this lncRNA at T5. The loss of iab6 lncRNA at the metasomal T5 may then interrupt the transcriptional or translational regulation of its target Hox genes, most likely Abd-B. Interruption of normal function of this Hox gene at T5 may then cause a homeotic shift of T5 setal identity toward the more anterior, black-colored T4 or T3, leading to the orange-to-black shift. Dashed arrows and question marks indicate untested alternative hypothesis on how ncRNA silencing may interrupt Abd-B during transcription and translation. EE, embryo epidermis; AE, adult epidermis.