Genome Architecture Facilitates Phenotypic Plasticity in the Honeybee (Apis mellifera)

Abstract

Phenotypic plasticity, the ability of an organism to alter its phenotype in response to an environmental cue, facilitates rapid adaptation to changing environments. Plastic changes in morphology and behavior are underpinned by widespread gene expression changes. However, it is unknown if, or how, genomes are structured to ensure these robust responses. Here, we use repression of honeybee worker ovaries as a model of plasticity. We show that the honeybee genome is structured with respect to plasticity; genes that respond to an environmental trigger are colocated in the honeybee genome in a series of gene clusters, many of which have been assembled in the last 80 My during the evolution of the Apidae. These clusters are marked by histone modifications that prefigure the gene expression changes that occur as the ovary activates, suggesting that these genomic regions are poised to respond plastically. That the linear sequence of the honeybee genome is organized to coordinate widespread gene expression changes in response to environmental influences and that the chromatin organization in these regions is prefigured to respond to these influences is perhaps unexpected and has implications for other examples of plasticity in physiology, evolution, and human disease.

Keywords: Apis mellifera; ChIP-seq; H3K27me3; RNA-seq; evolution of gene expression; gene cluster; genomic regulatory domain; honeybee; phenotypic plasticity; polycomb; reproductive constraint.

Fig. 1.

RNA-seq analysis of gene expression in the ovaries of queen-right workers (reproductively inactive), queen-less workers, and queens (both reproductively active). (A) Venn diagrams illustrate the number of genes differentially expressed in each pairwise comparison of the RNA-seq data (FDR corrected P value ≤ 0.01). (B) Gene ontology categories significantly enriched in queen-right ovaries (magenta) and queen-less ovaries (green), bars indicate number of genes, whereas spots indicate enrichment score. (C) Gene interaction network depicting physical and genetic interactions with the Notch signaling receptor in Drosophila melanogaster. Nodes are colored according to expression in honeybee ovaries (magenta, higher expression in queen-right; green, higher expression in queen-less; dark gray, not differentially expressed; light gray, not expressed in ovary; white, no known honeybee ortholog of the Drosophila gene). (D) Differential expression of polycomb group (PRC1 and PRC2) and trithorax group proteins (TAC1) in the honeybee ovary. Green, higher expression in queen-less (reproductively active) workers; gray, not differentially expressed.

Fig. 2.

Evolutionary histories of coordinately regulated clusters of differentially expressed genes within the Hymenoptera. (A) Heatmap depicting the conservation of “Queen responsive” clusters of coordinately regulated genes that are more highly expressed in queen-right (reproductively inactive) worker ovaries. (B) Heatmap depicting the conservation of “plasticity responsive” clusters of coordinately regulated genes that are more highly expressed in queen-less (reproductively active) worker ovaries. Intensity of the color depicts level of conservation along a continuous scale as indicated by the key. Clusters are ordered along the y-axis by the number of genes contained within the cluster the largest cluster the top of the heatmap. Phylogeny of the species is depicted on the y-axis (Peters et al. 2017). Social complexity is indicated by color of the species names (cyan, ancestrally solitary; green, facultative simple eusociality; orange, obligate simple eusociality; magenta, obligate complex eusociality). Species names are Amel, Apis mellifera; Emex, Eufriesea mexicana; Bimp, Bombus impatiens; Bter, Bombus terrestris; Mqua, Melipona quadrifasciata; Hlab, Habropoda laboriosa; Mrot, Megachile rotundata; Dnov, Dufourea novaeangliae; Hsal, Harpegnathos saltator; Pbar, Pogonomyrmex barbatus; Acep, Atta cephalotes; Aech, Acromyrmex echinatior; Pcan, Polistes canadensis; Nvit, Nasonia vitripennis.

Fig. 3.

Expression of genes of the PRC2 and localization of H3K27me3 in the honeybee ovary. (A) Ovary activity in queen-less worker bees is scored on a modified Hess scale. (B) RT-qPCR of genes encoding proteins of the PRC2. Target gene expression is measured relative to mRPL44 and Rpn2, which are stably expressed in honeybee ovaries (Duncan et al. 2016). Gene expression was measured in three biological replicates each consisting of ovaries from multiple individuals: queen (n = 3), queen-right worker (n = 20), and queen-less workers (score 0, n ∼ 20; score 1, n ∼ 20; score 2, n ∼ 10; and score 3, n ∼ 10). Differences in gene expression were assessed using a general linear model ANOVA with a Tukey post hoc test and 95% confidence interval. Samples that do not share letters are statistically significantly different with a P value <0.05. (C) Western blot analysis of ovary histone extracts for enrichment of H3K27me3 in queen, queen-right worker, and queen-less worker ovary (full blot: supplementary fig. 5, Supplementary Material online). (D) Inhibition of histone methylation using DZNep enhances ovary activity in honeybee workers. Proportion of bees scored as reproductively inactive (score, 0), and degrees of reproductively active (score 1–3) following treatment of newly emerged bees for 10 days with 50 μM DZNep (n = 524) or control (n = 532). Experiments were performed in triplicate on two separate occasions.

Fig. 4.

H3K27me3 stably defines “plasticity responsive” gene clusters, prefiguring changes in gene expression. (A) Average H3K27me3 enrichment across gene clusters more highly expressed in queen-right workers (Queen responsive clusters). Cyan is queen H3K27me3 enrichment, green is from queen-less workers, and magenta is from queen-right workers. (B) Boxplot illustrating H3K27me3 enrichment across gene clusters more highly expressed in queen-right workers (Queen responsive clusters). Only the 3′ flank of the clusters is significantly enriched for H3K27me3 in queen-right worker ovaries (Wilcoxon rank sum test, *P < 0.05, **P < 0.01, ns = not significant), showing that H3K27me3 marks are dynamic with respect to the presence of a queen and her pheromone. (C) Average H3K27me3 enrichment across gene clusters more highly expressed in queen-less worker ovaries (Plasticity responsive gene clusters) showing a decrease in H3K27me3 enrichment demarking both the 5′ and 3′ edges of the cluster. (D) Boxplot illustrating H3K27me3 enrichment across gene clusters more highly expressed in queen-less workers (Plasticity responsive clusters). In this case, we see that H3K27me3 marks are relatively stable, with both the 5′ and 3′ flanks of the cluster showing significant enrichment for H3K27me3 relative to the body of the cluster (Wilcoxon rank sum test, *P < 0.05, **P < 0.01, ns = not significant) even though the expression of these genes is low in queen-right workers. Boxplot whiskers indicate minimum and maximum, the box is defined by 25th percentile, median, and 75th percentile. Outliers, data points outside 1.5 times the interquartile range above the upper quartile and below the lower quartile, are shown as individual data points.

Fig. 5.

Network analysis identifies key hubs genes associated with differential enrichment of H3K27me. Genes with higher H3K27me3 enrichment in queen-right workers (magenta) (A) and queen-less workers (dark green) (B) are indicated by darker colored nodes. Interacting genes (queen-right workers, light magenta; queen-less workers, light green) were identified using BioGRID in DAVID. Network analysis was performed using Cytoscape. The predicted key hubs in this network have a high degree of centrality (as indicated by the relative size of each node) and are putative key regulators of reproductive constraint in honeybees.