DOI: 10.5281/zenodo.15428750
title: Intragenomic conflict associated with extreme phenotypic plasticity in queen-worker caste determination in honey bees (Apis mellifera)

description: Caste determination of honey bees (Apis mellifera) exemplifies developmental plasticity, where differences in larval diet result in identical genotypes yielding either long-lived, reproductive queens or short-lived, facultatively sterile workers. Beyond environmental factors, intragenomic conflict between genes inherited from the mother (matrigenes) versus the father (patrigenes) is hypothesized to influence underlying traits. In honey bees, the Kinship Theory of Intragenomic Conflict predicts selection on patrigenes to favor traits enhancing individual fitness, such as accelerated growth or increased body size – traits more active in queen-destined larvae. Although parent-of-origin effects have been reported in honey bees, the associated gene regulatory mechanisms remain unknown. While DNA methylation underlies parent-biased transcription in eutherians, angiosperms, and some insects, it is not operating in honey bees or other social insects. Allele-specific transcriptome analyses on queen- and worker-destined larvae at 192 hours post-fertilization (hpf) reveal hundreds of genes with parent-of-origin effects, with queen-destined larvae showing overrepresentation of patrigene-biased transcription relative to worker-destined larvae. Genes with parent-biased transcription resemble imprinted genes in other taxa regarding genomic clustering, recombination rate, intron length and CpG density, with a subset showing parent-biased transcription in 24 hpf eggs. Allele-specific ChIP-seq analyses demonstrate that parent-of-origin effects on caste-specific profiles of H3K27me3, H3K4me3 and H3K27ac are associated with parent-of-origin transcription. These findings suggest that parent-of-origin intragenomic conflict may influence phenotypic plasticity and may be associated with histone post-translational modifications, suggesting a “noncanonical” genomic imprinting-like system in social insects.
description: Sample CollectionReciprocally crossed colonies were generated from F1 male drones and female queens derived from different founder-generation (F0) stocks of European honey bees (Apis mellifera). From these crosses, we collected F2 192-hours post-fertilization (hpf) queen- and worker-destined larvae and 24 hpf eggs. F1 queens were produced by grafting, a standard commercial practice where young female larvae (72-96 hpf) are transferred to artificial queen cells containing royal jelly. These cells were placed in queenless colonies to induce queen-rearing behavior, then incubated at 33.5°C and 50% humidity until emergence. F1 drones were carefully managed to prevent colony drift, with newly emerged males marked with paint to identify their origin. Reciprocal crosses were created via instrumental insemination between 10-day-old F1 queens and 19-day-old F1 drones from different colonies within each block. Each inseminated queen was established in a mini colony of approximately 1,000 workers with queen excluders preventing further mating. Three weeks after queens began laying, labeled frames were provided for egg collection. From Blocks 1-2, some larvae were grafted at 96 hpf to induce queen development, with both queen- and worker-destined larvae collected at 192 hpf, frozen on dry ice, and stored at -80°C. From Blocks 3-4, pools of 24 hpf eggs (n=40 per pool) were collected and immediately preserved in TRIzol for later analysis.DNA Isolation & WGS Library Prepration & SequencingFlight muscle tissue was dissected from thoraces of F1 queen and drone parents (n=16) on a dry ice platform to maintain tissue integrity. Genomic DNA was isolated using DNeasy Blood & Tissue Kit (Qiagen). WGS libraries were prepared from 500 ng DNA using NEBNext Ultra II DNA Library Prep kits with NEBNext Multiplex Oligos (NEB), PCR-amplified for 3 cycles, and purified with AMPure XP Beads (Beckman Coulter).

RNA Isolation & mRNA Library Preparation & SequencingFor 192 hpf larvae (n=40) and 24 hpf eggs (n=12 pools of 40 eggs each), whole specimens were dissolved in 500 μl TRIzol, followed by chloroform extraction, isopropanol precipitation, and TURBO DNase (Invitrogen) treatment. RNA integrity was verified using TapeStation 4150 (Agilent), requiring RINe ≥9.0. For 192 hpf larvae, libraries were prepared from 1 μg total RNA using Illumina Stranded mRNA Prep kits with IDT for Illumina RNA UD Indexes (Set A), PCR-amplified for 10 cycles, and purified using AMPure XP Beads. For 24 hpf eggs, libraries were prepared from 500 ng RNA using the same kits, PCR-amplified for 10 cycles, and pooled at equimolar concentrations. All RNA libraries underwent quality assessment by D1000 DNA ScreenTape Analysis, dilution to ~2 nM, quantification using NEBNext Library Quant Kit on a QuantStudio 5 Real-Time PCR System, and treatment with Illumina Free Adapter Blocking Reagent before final cleanup.

Chromatin Preparation, ChIP & Library Preparation & SequencingChromatin was extracted from pooled 192 hpf F2 larvae (n=16 pools of 6 larvae). Samples were disaggregated in PBS with 5% Complete EDTA-free protease inhibitor cocktail (PIC) (Roche) using pellet pestles and needle passage. Cell viability was assessed with trypan blue staining. Approximately 30M cells per pool were cross-linked in 1% formaldehyde for 8 minutes, quenched with glycine, and lysed in a two-step process (cell lysis buffer followed by nuclear lysis buffer). Chromatin was sheared to 100-600 bp fragments (average 250 bp) using a Bioruptor Pico (Diagenode) for 5 cycles. For immunoprecipitation, antibody-bead complexes were prepared using 5 μl (1 μg/μL) of antibodies against H3K27me3 (Abcam #ab6002), H3K27ac (Active Motif #39133), or H3K4me3 (Abcam #ab12209) with 25 μl Protein G Dynabeads (Invitrogen). Beads were blocked with BSA, combined with 450 μl (~12 μg) pre-cleared chromatin, and incubated for 6 hours. After sequential washes with buffers B, C, D, and TE, bound chromatin was eluted and processed to remove RNA contamination and reverse cross-links. ChIP-seq libraries (n=56 total, including 8 input controls) were prepared using NEBNext Ultra II DNA Library Prep kits. Input libraries were prepared from 200 ng DNA with 5 PCR cycles, while ChIP libraries used 1-10 ng DNA with 12 PCR cycles. All libraries underwent cleanup, quality assessment, quantification, and pooling procedures as described for RNA libraries.Summary of Statistical AnalysesFor differential gene expression analysis, we used Kallisto to estimate transcript abundances from larval samples, with transcripts combined into a gene expression matrix using tximport. Sample relationships were assessed via k-means clustering, and DESeq2 was employed to test for differential expression between castes (queen vs. worker), treating lineage as a cofactor. Effect sizes were adjusted using Bayesian shrinkage, with focus on differentially expressed genes at thresholds of |log2(fold change)| > 1.5 and FDR < 1%.To identify single nucleotide polymorphisms (SNPs), whole-genome sequencing reads were trimmed, aligned to the A. mellifera reference genome, and filtered to remove low-quality alignments. Variants were detected using freebayes and filtered for quality. High-confidence homozygous SNPs were integrated into the reference genome and intersected to identify parent-specific SNPs, which were then mapped to gene transcripts.For allele-specific transcriptome analysis, mRNA-seq reads were aligned to parent genomes, and coverage was calculated at each SNP. We employed both the Storer-Kim binomial exact test and a general linear mixed effects model to identify transcripts with parent-biased expression. To control for sequencing depth variation, we normalized read counts and used size factors as offset terms. Significant transcripts had to show consistent bias direction across SNPs with specific thresholds for maternal bias (p1 < 0.4 & p2 > 0.6), paternal bias (p1 > 0.6 & p2 < 0.4), lineage B bias (p1 < 0.4 & p2 < 0.4), or lineage A bias (p1 > 0.6 & p2 > 0.6). Chi-squared tests determined whether queen-destined larvae showed enrichment for specific allelic biases compared to worker-destined larvae.Weighted Gene Co-expression Network Analysis (WGCNA) constructed unsigned co-expression network modules using the variance-stabilized gene expression matrix. Modules with eigengene correlation coefficients |r| > 0.75 were merged, and correlation with caste was calculated. Hypergeometric tests identified modules enriched for genes with parent-biased transcription.For ChIP-seq analysis of histone modifications, reads were trimmed, aligned to the reference genome, and filtered for quality. Problematic regions were excluded using a blacklist of repetitive and high-noise regions. Broad and narrow peaks were called using MACS2, and consensus peaks were identified for each histone modification by caste. Reads were aligned to parent genomes to identify allelic enrichment within consensus peaks, with separate analyses for promoters and gene bodies.Genome-wide enrichment differences for histone modifications were assessed using csaw for sliding windows across the genome. Statistical tests of differential enrichment between castes used edgeR, treating lineage as a cofactor.Logistic regressions evaluated relationships between histone modification enrichment and transcription, as well as between parent-biased transcription and parent-biased histone modification enrichment.Gene ontology enrichment analyses were conducted with DAVID, focusing on biological processes and KEGG pathways at FDR < 5%. Sequence feature analyses examined spatial clustering of genes, intron length, GC content, CpG density, and recombination rates, comparing parent-biased and unbiased genes using independent t-tests.