Beyond conflict: Kinship theory of intragenomic conflict predicts individual variation in altruistic behaviour

Behavioural variation is essential for animals to adapt to different social and environmental conditions. The Kinship Theory of Intragenomic Conflict (KTIC) predicts that parent‐specific alleles can support different behavioural strategies to maximize allele fitness. Previous studies, including in honey bees (Apis mellifera), supported predictions of the KTIC for parent‐specific alleles to promote selfish behaviour. Here, we test the KTIC prediction that for altruism‐promoting genes (i.e. those that promote behaviours that support the reproductive fitness of kin), the allele with the higher altruism optimum should be selected to be expressed while the other is silenced. In honey bee colonies, workers act altruistically when tending to the queen by performing a ‘retinue’ behaviour, distributing the queen's mandibular pheromone (QMP) throughout the hive. Workers exposed to QMP do not activate their ovaries, ensuring they care for the queen's brood instead of competing to lay unfertilized eggs. Due to the haplodiploid genetics of honey bees, the KTIC predicts that response to QMP is favoured by the maternal genome. We report evidence for parent‐of‐origin effects on the retinue response behaviour, ovarian development and gene expression in brains of worker honey bees exposed to QMP, consistent with the KTIC. Additionally, we show enrichment for genes with parent‐of‐origin expression bias within gene regulatory networks associated with variation in bees' response to QMP. Our study demonstrates that intragenomic conflict can shape diverse social behaviours and influence expression patterns of single genes as well as gene networks.

same stimuli (Laskowski et al., 2022).One of the most complex behaviours is communication, which requires the production, reception and interpretation of signals which in turn influence the behaviour of individuals (Wilson, 1975).In social groups, variation in receptivity to these signals underlies the flexible task repertoires which facilitate collective responses to changing environments (Gordon, 2014(Gordon, , 2016)).
Studying the proximate and ultimate mechanisms underlying variation in complex behaviours such as signal receptivity among social group members can shed light on the adaptive significance of these traits, how they arise and why and how they persist in different populations.
The Kinship Theory of Intragenomic Conflict (KTIC) provides a fascinating framework in which to examine how proximate and ultimate mechanisms interact to influence gene expression and phenotypic outcomes (Gardner & Úbeda, 2017;Haig, 1992Haig, , 2000)).The KTIC proposes that alleles derived from the mother (matrigenes) versus the father (patrigenes) can have different strategies for maximizing their reproductive fitness-with conflicting effects on development and behaviour-and these different selective pressures lead to opposite expression patterns of matrigenes and patrigenes.Since these selective advantages and expression differences arise from the origin of the gene (whether it is derived from the mother or the father) and not the gene sequence, understanding how these parentof-origin expression differences are regulated and selected for represents a fascinating question in molecular ecology (Ferguson-Smith, 2011;Oldroyd & Yagound, 2021).
Intragenomic conflict related to social behaviour is predicted to occur frequently in species in which offspring are asymmetrically related-when matrigenes and patrigenes are not distributed in the same way (Bourke, 2014;Gardner & Úbeda, 2017;Haig, 1992Haig, , 2000)).If females are polyandrous, their offspring share matrigenes but have different patrigenes, and thus, there is a relatedness asymmetry between offspring of different fathers.In such organisms, when offspring use maternal resources for survival, matrigenes should be selected to act altruistically by reducing the use of maternal resources and increasing cooperation between offspring.In this context, altruism ensures that sufficient resources are available for future offspring, which have a higher probability of sharing matrigenes than patrigenes.In contrast, patrigenes are predicted to act selfishly by increasing the use of maternal resources and competition between offspring, which have a lower probability of sharing patrigenes than matrigenes (Patten et al., 2014).Previous empirical tests of the KTIC have focussed on genes in which patrigenes favour acting 'selfishly', such as during mammalian embryonic development and fetal growth (Barlow & Bartolomei, 2014;Moore et al., 2015), seed production in plants (Rodrigues & Zilberman, 2015) and reproduction in worker honey bees (Galbraith et al., 2016(Galbraith et al., , 2021)).The KTIC also predicts that for altruism-promoting genes, the allele with the higher altruism optimum (the matrigene in the examples above) should be selected to be expressed while the other is silenced (Hamilton, 1964), but empirical tests of this hypothesis are still lacking.
Honey bees are an outstanding model system for studying the causes and consequences of intragenomic conflict (Bresnahan et al., 2023;Galbraith et al., 2016Galbraith et al., , 2021;;Gibson et al., 2015;Kocher et al., 2015;Smith et al., 2020).Due to their haplodiploid genetics, polyandrous mating system and large colony size, matrigenes and patrigenes are unevenly distributed across colony members, and multiple social interactions are predicted to be influenced by intragenomic conflict (Queller, 2003).Honey bee colonies have a single reproductive female queen, while her daughters (workers) typically remain sterile and forego their own reproduction to rear the queen's offspring (Page, 2013).If the queen is lost, workers can act selfishly, activating their ovaries and competing among each other to lay unfertilized eggs, or altruistically, by remaining sterile and rearing a new queen (Galbraith et al., 2015).Since queens mate with multiple males, daughters on average share more of their matrigenes than patrigenes: Thus, patrigenes benefit from workers activating their own ovaries and laying eggs while matrigenes benefit from workers remaining sterile and helping rear the queen or their sister's offspring (Queller, 2003).Indeed, previous studies in honey bees demonstrated that patrigenes favour the selfish behaviours of worker ovary activation (Galbraith et al., 2016(Galbraith et al., , 2021) ) and worker ovary development (Smith et al., 2020).
Honey bees are also an excellent system in which to evaluate whether and how intragenomic conflict can lead to variation in altruistic behaviour, including in their receptivity to a social signal.
Worker bees act altruistically when they forgo their own personal reproductive fitness in favour of queen and sibling brood care (Ratnieks & Helanterä, 2009).A pheromone produced by the queen-in particular, a five-component blend termed queen mandibular pheromone (QMP)-triggers behavioural and physiological changes in workers (Grozinger, 2015;Traynor et al., 2014).Bees exposed to QMP show reduced worker ovary activation (Mumoki et al., 2018) and extended brood care activities (Pankiw et al., 1998).Moreover, QMP is a contact pheromone that releases a 'retinue response' behaviour in workers, attracting them to surround and groom the queen as she moves across the brood comb (Keeling et al., 2003).
While tending to the queen in the retinue response, worker bees pick up QMP and spread it to other workers throughout the colony (Seeley, 1979).Thus, the retinue response is an altruistic behaviour, whereby the workers increase their exposure to QMP and spread it to the other bees, ensuring that all workers remain sterile and rear the queen's offspring.According to the KTIC, due to the haplodiploid genetics of honey bees and multiple mating by the queen, the retinue response behaviour should be favoured by matrigenes (Queller, 2003).While several studies have shown that there can be significant variation in retinue response behaviour (Galbraith et al., 2015;Kocher et al., 2010), the role of intragenomic conflict has not been investigated.
In this study, we tested predictions of the KTIC on the altruistic, pheromone-mediated retinue response behaviour in worker bees.
Using a reciprocal cross design, we assessed parent-of-origin effects on, and relationships between, responsiveness to QMP, ovary activation and transcriptional profiles in brains of worker bees.First, we tested the hypothesis that the physiological effects of QMP on workers (i.e.smaller, and less active ovaries) are influenced by a worker's parent-of-origin.Second, we conducted an allele-specific transcriptome analysis of worker bees that differed in QMP responsiveness to test the hypothesis that a worker's response to QMP is associated with an enrichment for maternally biased gene expression.Specifically, we hypothesized that workers that responded to QMP by performing the retinue behaviour should show more genes with maternal allele-biased expression compared with workers that did not respond to QMP.Although QMP may exert some of its effects by causing workers to silence certain genes-for example, those involved in ovary development-a set of proteins must be required to mediate the reception of and response to QMP; therefore, we expected more of such genes to show maternally biased expression in workers that responded to QMP compared with those that did not.
These studies also allowed us to examine the question of how parent-specific gene expression leads to individual variation in behaviour and fitness.The KTIC predicts that selection should favour a balance between the individual effects of matrigenes and patrigenes engaged in conflict (Haig, 2008;Patten et al., 2016).In support of this prediction, previous studies in honey bees-which demonstrated that patrigenes favour the selfish behaviours of worker ovary activation (Galbraith et al., 2016(Galbraith et al., , 2021) ) and aggression (Bresnahan et al., 2023)-found that genes showing differences in parent-specific transcript abundance do not show differences in overall transcript abundance, presumably because upregulation of one allele was balanced by downregulation of the other allele.If the total amount of transcript for a given gene is the same between two behavioural phenotypes, how do genes engaged in parent-of-origin conflicts contribute to variation in physiology or behaviour?It has been hypothesized that genes engaged in parent-of-origin conflicts exert their influence on behavioural phenotypes through gene networks (Patten et al., 2016).In support of this hypothesis, genes showing intragenomic conflict are frequently enriched in gene networks (Al Adhami et al., 2015;Bonthuis et al., 2022;Galbraith et al., 2021).Additionally, allelic differences in the occupancy of transcription factors (TFs)-the regulatory component of gene networks-are associated with allelic differences in transcription (Do Kim et al., 2006;McDaniell et al., 2010;Reddy et al., 2012).Therefore, small changes in expression of matrigenes and patrigenes, potentially mediated by differences in transcription factor (TF) activity, may result in downstream effects on gene networks.Here, we explored how genes that show intragenomic conflict are associated with gene regulatory networks (GRNs) underlying behavioural variation.Specifically, we hypothesized that gene networks that are correlated with-and genes that are predictive of-the behavioural state of worker bees exposed to QMP, should be enriched for genes, and the targets of TFs, which exhibit parent-of-origin expression bias.

| Biological samples
Reciprocally crossed colonies were generated between F0 stocks of Apis mellifera ligustica × scutellata hybrid honey bees (i.e.'Africanized' honey bees, or 'AHBs') and European honey bees ('EHBs'; derived from A. mellifera subspecies that evolved in Europe) and between different stocks of EHBs.In studies conducted in 2019, we sourced reproductive females (queens) and males (drones) from two AHB colonies (Colonies

F I G U R E 1 Reciprocal cross design.
From each source colony (C1-C12), queens (♀) and drones (♂) were selected and crossed (AxB, BxA) to construct reciprocal crosses (RC1-RC14).Reciprocal crosses that were selected for sequencing are underlined.Colours indicate bee stock.Blocks 1-2 were generated for studies conducted in 2019, Blocks 3-4 for studies in 2021, and Blocks 5-6 for studies in 2022.
From each F0 colony, F1 queens were generated using a standard commercial practice known as grafting (Connor, 2009).Adult virgin queens were labelled with a numbered tag affixed to the thorax and housed in cages within a nursery colony until F1 drones reared from these colonies reached sexual maturity.As drones tend to drift between colonies (Currie & Jay, 1991), frames of emerging F1 drone brood were identified and isolated between queen excluders to prevent escape.Newly emerged F1 drones were marked on the thorax with nontoxic paint to identify their colony of origin and released back into the main colony entrance.In Blocks 1-4, two F1 queens and two F1 drones were selected from each colony and crossed with two F1 queens and two F1 drones from the reciprocal colony by instrumental insemination (Cobey et al., 2013) to create two reciprocal crosses per block (reciprocal crosses 1-8, or RC1-RC8, Figure 1).In Blocks 5 and 6, the same procedure was conducted with three F1 queens and three F1 drones from each colony to create three reciprocal crosses per block (RC9-RC14, Figure 1).Inseminated queens were placed in their own colonies to initiate egg-laying, and the entrances of these colonies were restricted with queen excluder material to prevent escape or further mating.
Approximately 3 weeks after the F1 queens began laying fertilized eggs, frames of sealed, emerging F2 brood were collected from each colony and stored in separate containers within a standing incubator kept at constant temperature of 34.5°C and 75% relative humidity.
Once collection of the emerging F2 brood was completed, the queen from each colony was euthanized and stored at -80°C for DNA extraction and sequencing.The bodies of the drones used for inseminating each queen were likewise placed on dry ice and stored at -80°C immediately after semen collection for subsequent analysis.
In Blocks 1-4, approximately five hundred age-matched newly emerged F2 worker bees from each F1 cross were collected on Day 1 of the experiment and marked on the thorax with nontoxic paint to identify their maternal lineage of origin (Figure S1a).These bees were combined with an equal number of worker bees from the respective reciprocal crosses and placed into 'nucleus' hive boxes containing honeycomb frames with honey and pollen but no queen or brood.In total, we created eight small nucleus colonies, each containing 500 F2 worker bees from each respective F1 cross in a reciprocal cross design.
We also added 1000 unmarked worker bees from unrelated colonies of European origin to each hive to increase the number of worker bees in the nucleus colony.Worker bees were kept in a 'queenright' state (emulating exposure to the queen) by providing the colony with a plastic TempQueen strip (Mann Lake, Wilkes-Barre, PA) infused with QMP, as done previously (Pankiw et al., 1994).
In Blocks 5 and 6, for each reciprocal cross, we populated six cages with five number-tagged age-matched newly emerged F2 worker bees from each F1 cross (Figure S1b).Workers were housed in petri dishes with pieces of beeswax foundation (Brushy Mountain Bee Farm, Moravian Falls, NC) fit against the back wall, following methods described in Shpigler and Robinson (2015).Each cage was provided with a pollen ball at the bottom of the dish and a 50% sucrose solution via a modified 1.5 mL microcentrifuge tube for the bees to feed ad libitum.Worker bees were kept in a 'queenright' state by providing a removable glass slide dosed with 0.1 queen equivalents of fresh QMP (Intko Supply LTD, Vancouver, BC) diluted in 1% water/isopropanol changed daily, following methods described in Kocher et al. (2010).

| Parent-of-origin effects on the retinue response behaviour
In Blocks 1-4, we measured the colony-level frequency of retinue behaviour as the proportion of individuals from each cross that was observed performing the retinue response (i.e.touching with antennae and mouthparts) to the pheromone lure.On Day 6 of the experiment, we affixed a fresh TempQueen strip to the centre of an empty honeycomb frame (containing no food or brood) placed in the middle of each colony.Approximately 20 min after introducing the pheromone lure, we retrieved the frame.Over a 10-min period, we collected paint-marked worker bees observed antennating the pheromone lure as being in a 'responsive' behavioural state, and worker bees observed on the edges of the frame as being in an 'unresponsive' behavioural state (Table S1).To test for an effect of maternal lineage (source colony) on the proportion of worker bees from each reciprocal cross colony that were observed performing the retinue behaviour, we conducted logistic regression in R, treating the maternal lineage as a fixed effect nested within the block.We conducted this assay once for Blocks 1 and 2 and three times on the same day spaced 3 h apart for Blocks 3 and 4 (Table S2).
In Blocks 5 and 6, we measured the individual-level frequency of retinue response behaviour as in Kocher et al. (2010).On Days 4-6, we provided fresh QMP to each mini colony.Approximately 20 min after introducing the pheromone, we recorded responses to the stimulus (i.e.antennating and touching the slide where the QMP was added with their mouthparts) every 5 min for a total of 25 min, scoring a response as '1' and no response as '0' (Table S3).We used a linear mixed-effects model, treating the cage as a random effect, to test for an effect of maternal lineage nested within the block on the retinue response scores of individual worker bees.Additionally, we performed one-way ANOVA tests to determine whether the frequency of retinue response across days varied between cages, between individual bees within cages and/or between individual bees overall (Table S4).

| Parent-of-origin effects on ovary size and activation stage
The ovaries of 10 individual worker bees per behavioural state, per cross, per reciprocal cross and from the same blocks used for sequencing (10 bees × 2 behavioural states × 2 crosses × 2 reciprocally crossed colonies × 3 blocks = 240 individuals in total) were dissected under 1× phosphate-buffered saline.Dissected ovaries were placed on a microscope slide and observed at 5× magnification for measurements of ovary size (the number of ovarioles present) and activation stage (four discrete stages or scores based on a modified Hess scale (Hess, 1942)) of 0-3 as done previously (Duncan et al., 2016; Table S5).For each block, separately, one-way ANOVA tests were used to determine the significance of variation in ovary size and activation stage, respectively, by behavioural state and maternal lineage nested within behavioural state (Table S6).Additionally, one-way ANOVA tests were used to determine the significance of overall variation in ovary size and activation stage, respectively, by behavioural state and maternal lineage nested within block (Table S7).

| Sample preparation for whole genome resequencing and RNA-seq
The F1 queen and drone parents from each of three reciprocal crosses (RC2, RC5 and RC8, underlined in Figure 1) constituting three distinct genetic blocks (Blocks 1, 3 and 4) were selected for whole genome resequencing (WGS).The thorax of each bee was removed, and the flight muscle tissue was dissected on a platform surrounded by dry ice to keep the tissue frozen until the DNA isolation procedure.Genomic DNA was isolated using a DNeasy Blood & Tissue Kit (Qiagen, Germantown, MD).Samples were sent to the Penn State Genomics Core Facility for quantitation by a Qubit Fluorometer (Life Technologies, Carlsbad, CA).Sequencing libraries were prepared using an Illumina DNA prep kit with unique 5′ and 3′ indexes, pooled and sequenced on a NextSeq 2000 P3 flow cell to generate approximately 37 million 150-bp paired-end reads per sample (Table S8).S8).WGS and RNA-seq reads for this study are deposited in the NCBI Sequence Read Archive under SRA project accession PRJNA732718.

| Construction of F1 genomes, identification of SNPs within transcripts and detection of allele-specific gene expression in F2 workers
We utilized the Storer-Kim (SK) binomial exact test of two proportions (Storer & Kim, 1990;Wang & Clark, 2014) and a general linear mixed-effects model with interaction terms (GLIMMIX) to identify transcripts that exhibited parent-of-origin expression bias following methods described previously (Bresnahan et al., 2023;Galbraith et al., 2016Galbraith et al., , 2021;;Kocher et al., 2015).Briefly, WGS reads were aligned to the latest A. mellifera genome assembly (Wallberg et al., 2019) to identify SNPs, which were used to construct individual genomes for each parent.Offspring RNA-seq reads were then aligned to their respective parental genomes for counting parentof-origin reads within transcripts.SK and GLIMMIX tests were then conducted in R (R Core Team, 2022) using the allele-specific read counts.Full details of these methods are included in the Appendix S1 document associated with this manuscript.Supplementary datasets associated with this analysis are provided in Tables S9-S19.

| Differential gene expression analysis
We used Kallisto (Bray et al., 2016) to estimate transcript abundances.A transcriptome index was generated with gffread (Pertea & Pertea, 2020) using the Amel_HAv3.1 RefSeq genome feature file.Transcript abundance files were combined to generate a gene expression matrix (Table S18) using the tximport package (Soneson et al., 2022) in R for analysis of differential expression with DESeq2 (Love et al., 2014).After identifying and removing an outlier library for one individual bee using a PCA projection pursuit method (Chen et al., 2020; Figure S6), we employed a two-factor design to test for differential expression between behavioural states, treating maternal lineage as a cofactor (Table S20).

| Support vector machine (SVM) classification
Support vector classification (SVC), also known as SVM, is a machine learning tool used for binary classification tasks.The goal of an SVM is to identify a hyperplane that best separates data points belonging to different classes while maximizing the distance between the hyperplane and the nearest data points of each class.
Applied to transcriptomics data, SVMs have been used to classify cancer subtypes from tumour samples (Huang et al., 2018) and behavioural phenotypes of social insects (Taylor et al., 2021).
Here, SVM classification was conducted using the transcriptomes of QMP-unresponsive bees (coded with a value of 0) and QMPresponsive bees (coded with a value of 1).Transcriptomes for SVM classification were normalized to adjust for variation in sequencing depth between libraries by variance stabilizing transformation followed by quantile normalization.Genes with low counts and zero variance between samples were filtered from the dataset, leaving 11,514 of the 12,123 genes (94.98%) for SVM training.We subset the transcriptomic data on a random selection of samples into 80% training and 20% testing sets.Classifiers were fit using a radial kernel function, and the hyperparameters cost (C) and gamma (γ) were tuned via grid search in a parameter space of 2 −2 < C < 1 and 10 −3 < γ < 10 −1 .Feature selection was performed following methods described in Taylor et al. (2021).Starting with a classifier fit with all 11,514 genes, we iteratively performed the following: (1) the mean threefold cross-validation error of 20 rounds using random bins was calculated; (2) the matrix product of the classifier's coefficients and its support vectors was used to calculate the feature weights of all remaining genes; (3) the gene with the smallest absolute weight was removed, until only 100 genes (an arbitrary minimum number of features chosen a priori) remained.
The support vector classifier with the minimum cross-validation error or minimum number of features was then selected for downstream analyses (Figure S7; Table S21).

| Gene regulatory network inference
Weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA package (Langfelder & Horvath, 2008) in R. Signed co-expression network modules were constructed using a soft-threshold power of 15 to balance scale independence and mean connectivity (Figure S8).Modules with eigengene correlation coefficients >0.75 were merged (Figure S9; Table S22).We then calculated correlation coefficients of each module with behavioural state.Hypergeometric tests were conducted to identify modules enriched for genes with parent-of-origin expression bias and genes in our minimal support vector classifier of behavioural state.
The Analyzing Subsets of Transcriptional Regulators Influencing eXpression (ASTRIX) method for gene regulatory network (GRN) inference (Chandrasekaran et al., 2011;Jones et al., 2020;Shpigler et al., 2017Shpigler et al., , 2019) ) was performed in MATLAB.Lists of honey bee orthologs of Drosophila melanogaster TF genes were retrieved from Kapheim et al. (2015) and Jones et al. (2020).TFs predicted to interact with a target gene by ASTRIX were required to fulfil three criteria: (1) share a significant degree of mutual information with the target gene (p < 1 × 10 −6 ), (2) explain at least 10% of the variance of the target gene, quantified by least angle regression and (3) be predicted with a correlation of at least 0.8 using expression levels of the TF and target gene (Table S23).We then subset the GRN to include only those genes predicted to interact with parent-of-origin biased TFs (Table S24).With this subset of the GRN, we performed hypergeometric tests to identify gene co-expression network modules enriched for the targets of parent-of-origin biased TFs.

| Gene ontology (GO) enrichment analyses
GO enrichment analyses were conducted with DAVID (Sherman et al., 2022) using D. melanogaster orthologs of the genes showing parent-of-origin expression bias (Table S25), the genes in our minimal support vector classifier of behavioural state (Table S26) and the gene co-expression network modules constructed with WGCNA (Table S27).Drosophila melanogaster orthologs of A. mellifera genes were retrieved from Bresnahan et al. (2022).For this study, we focussed on GO terms (specifically, biological processes and molecular functions) and KEGG pathways found to be significantly overrepresented among candidate genes at a false discovery rate of < .05.

| Parent-of-origin effects on the retinue behaviour
Our results suggest a parent-of-origin effect on-and considerable individual variation in-the retinue response behaviour in worker bees.We observed a significant association between maternal lineage and the proportion of worker bees performing the retinue response behaviour in Blocks 3 and 4, but not in Blocks 1 and 2 (Figure 2a; Table S2).Additionally, in Blocks 5 and 6, we observed a significant association between maternal lineage and the frequency of retinue response across days of individual bees (Figure 2b), between individual bees (Figure 2c), between cages (Figure S2) and between individual bees within cages (Table S4).

| Parent-of-origin effects on ovary activation
To test for an association between maternal lineage and reproductive physiology, we evaluated ovary size and activation in Blocks 1, 3 and 4, which included the reciprocal crosses that we subsequently used for sequencing.Bees in Blocks 5 and 6 were not used for the remainder of our study as they were collected as older nurses for a separate study.We found that workers from colonies showing more frequent responses to QMP had smaller ovaries in all tested blocks, but only those in Blocks 3 and 4 showed a significant effect of maternal lineage.Specifically, individuals from the A ♀ B ♂ cross in Blocks 3 and 4 had larger ovaries (Figure 3a, Table S6)-notably, these individuals also responded less frequently to QMP (Figure 2a).In a full model including all three blocks, ovary size varied significantly by behavioural state and maternal lineage (Table S7).These data suggest there is a parent-of-origin effect on ovary size and there are larger ovaries in QMP-unresponsive compared with QMP-responsive bees.
Ovary activation varied significantly by behavioural state in all tested blocks and by maternal lineage in Blocks 3 and 4 but not in Block 1 (Figure 3b, Table S6).Specifically, in all three blocks, ovaries were more active in QMP-unresponsive compared with QMPresponsive bees, and ovaries were more active in individuals from the A ♀ B ♂ cross in Blocks 3 and 4. In a full model including all three blocks, ovary activation varied significantly by behavioural state and maternal lineage (Table S7).These data suggest there is a parent-oforigin effect on ovary activation, and bees that did not respond to QMP had more active ovaries.Taken together, we found that bees from the maternal lineage that responded less frequently to QMP had larger and more active ovaries.

| Parent-of-origin gene expression
At approximately 40× genome coverage, we detected an average of 2.77 million homozygous SNPs per F1 diploid queen and 2.87 million SNPs per F1 haploid drone.On average, 550,695 of these SNPs were within transcripts and varied in their sequence identity between the parents of each cross, allowing for the identification of parent-of-origin reads in the F2 worker bees, 491,338 (89%) of which were within transcripts that had at least two SNPs.After filtering SNP positions with low RNA-seq read coverage in the F2 bees, our datasets contained read coverage at an average of 77,131 positions distributed among an average of 4366 of 11,865 transcripts (36.8%) in 12 QMP-unresponsive bees per block, and an average of 78,063 positions distributed among an average of 4395 of 11,865 transcripts (37%) in 12 QMP-responsive bees per block (Tables S9-S11).Considering the overlap between behavioural states, on average, bees in each block expressed 4143 of 11,865 transcripts (34.92%).
Accounting for all unique transcripts across the three blocks, we identified 624 maternally biased transcripts in whole brains of unresponsive bees and 887 in responsive bees, in addition to 143 paternally biased transcripts in unresponsive bees and 181 in responsive bees (Table S19).Of these transcripts, 64 showed consistent allele-biased expression in at least two blocks, all of which were biased towards the  S25).

| Gene expression differences
In both a full model accounting for all three blocks and a reduced model excluding the blocking factor, we did not detect any significantly differentially expressed genes in whole brains of QMPunresponsive worker bees compared with QMP-responsive bees (Table S20).We therefore used SVC to identify more subtle gene expression differences which best classify individuals into behavioural states (Taylor et al., 2021).We trained a SVM using transcriptomes of 14 QMP-responsive bees and 15 QMP-unresponsive bees.A full model based on all 11,514 genes which passed our low count threshold achieved a root mean squared validation error of 0.5105 in threefold cross-validation.To identify model fit with a minimal set of genes that were maximally predictive of behavioural state, we conducted feature selection to progressively remove uninformative genes.The minimal model, which was trained on 100 genes (Table S21), had a root mean squared classification error of 0.2093, which was a significant improvement over that achieved by the full model (Figure S5).Genes in the minimal model were enriched for 15 GO terms, including negative regulation of histone modifications, negative regulation of DNA recombination, nucleosome assembly and chromosome organization (Table S26).Parentof-origin biased genes were not significantly overrepresented in the model (8/100, p = .94).

| Inference of gene regulatory networks
We performed weighted gene co-expression network analysis (WGCNA) to identify modules of co-expressed genes (Langfelder & Horvath, 2008).Additionally, we used the ASTRIX method (Chandrasekaran et al., 2011;Jones et al., 2020;Shpigler et al., 2017Shpigler et al., , 2019) ) to determine whether intragenomic conflicts within these coexpressed gene modules are influenced by TFs with parent-of-origin expression bias.For both analyses, we used the same 11,514 genes adjusted by variance stabilizing transformation and quantile normalization as in SVC.

F I G U R E 4
Worker bee response to QMP is associated with maternal allele-biased transcription.Allele-specific transcriptomes were assessed in worker bees from a reciprocal cross between different stocks of EHBs (C5 × C9, Block 3) that were unresponsive or responsive to QMP.The x-axis represents, for each transcript, the proportion of cross A (C9 lineage) reads in bees with a cross B (C5 lineage) mother and cross A father (p1).The y-axis represents, for each transcript, the proportion of cross A reads in bees with a cross A mother and cross B father (p2).Each colour represents a transcript which is significantly biased at all tested SNP positions: black is maternal, purple is cross A, gold is cross B, blue is paternal, and grey is not significant.Significance was determined using the overlap between two statistical tests: a generalized linear interactive mixed model (GLIMMIX), and a Storer-Kim test along with previously established cut-off thresholds of p1 < .4 and p2 > .6 for maternal bias, p1 > .6 and p2 < .4 for paternal bias, p1 < .4 and p2 < .4 for lineage B bias and p1 > .6 and p2 > .6 for lineage A bias.Similar results were obtained in Blocks 1 and 4, see Figures S3 and S4 for details.
Using this approach, we identified 22 co-expressed gene modules (Table S22) and inferred a GRN consisting of 4200 interactions between 196 TFs and 1839 target genes (Table S23).Within this network, 29 of the 196 TFs (14.8%) showed parent-of-origin expression bias and were predicted to interact with 353 of the 1839 target genes (19.2%).Specific co-expression modules were (1) correlated with behavioural state, (2) overrepresented for parent-biased genes, (3) overrepresented for genes in our minimal support vector classifier of behavioural state and (4) overrepresented for the targets of parent-biased TFs (Table 1).
Two co-expression modules (Modules 4 and 7) were significantly correlated with behavioural state.Module 4 was enriched for genes in the Notch signalling pathway (Table S27), which has been demonstrated previously to mediate worker reproductive constraint in the presence of the queen (Duncan et al., 2016).Module 7 was not enriched for any GO terms or KEGG pathways, but it was overrepresented for pheromone binding proteins and signalling peptides with known roles in honey bee behavioural maturation, including vitellogenin (Vg), juvenile hormone esterase (JHe) and hexamerin 70a.
Six co-expression modules were significantly overrepresented for genes with parent-of-origin expression bias, including Modules 4 and 7 described above.Three of these additional modules (Modules 2, 3 and 8) were enriched for KEGG pathways and/or GO terms (Table S27).Module 2 was enriched for genes in the MAPK and Hippo signalling pathways, dorso-ventral axis formation, in addition to 109 GO terms for neural development, sensory perception, circadian rhythm, transcription regulation and chromatin organization.Module 3 was enriched for genes in the Wnt signalling pathway, spliceosome, neuroactive ligand-receptor interactions, in addition to 76 GO terms for transcription regulation, ion transport, synaptic transmission and neural development.Additionally, Module 8 was enriched for genes with protein binding function.Finally, three co-expression modules were enriched for genes predicted to be regulated by parent-biased TFs (Figure S6c), including Modules 2 and 3 described above, in addition to Module 10, which was enriched for genes related to nucleosome assembly and DNA-templated transcription (Table S27).
Notably, Module 10 was also enriched for genes in our minimal support vector classifier of behavioural state (Table 1).

TA B L E 1
Gene co-expression network modules are significantly correlated with behavioural state and enriched with parent-biased genes, behaviour-state classifying genes and targets of parent-biased TFs.

| Variation in altruistic behaviour is associated with variation in matrigene expression
Here, we report evidence of a parent-of-origin effect on an altruistic behaviour and associated physiological traits and on transcription profiles of whole insect brains.Specifically, we found the proportion of worker bees within a colony responding to QMP was associated with the bees' maternal lineage, as was the frequency of retinue response between individual bees (Figure 2; Figure S2; Table S4).
Ovary size and activation were also associated with maternal lineage, with bees that had larger and more active ovaries being less responsive to QMP (Figures 2 and 3; Tables S6 and S7).In each reciprocal cross colony, bees had larger and more active ovaries when they came from crosses that were less responsive to QMP.Moreover, we found that maternally biased transcripts were enriched in QMP-responsive worker bees relative to QMP-unresponsive bees (Figure 4; Figures S3 and S4).These results support a model in which bees with smaller ovaries-a trait which is determined during preadult development (Hartfelder et al., 2018)-are less likely to successfully compete with their sisters for ovary activation if the queen is lost (Makert et al., 2006) and are more likely to engage in the retinue response, increasing both their own and their sisters' exposure to queen pheromone (Page, 2013;Seeley, 1979).Together with previous studies, these results indicate that matrigenes favour the development of smaller ovaries while patrigenes favour the development of larger ovaries during preadult development.Likewise, matrigenes favour the retinue response while patrigenes inhibit it, and matrigenes inhibit ovary activation while patrigenes favour it (Galbraith et al., 2015(Galbraith et al., , 2016;;Kocher et al., 2015;Kocher & Grozinger, 2011;Reid et al., 2017;Smith et al., 2020).Thus, during an individual worker's life, patrigenes and matrigenes are in conflict over developmental, physiological and behavioural traits.
One point of contention in our results is the absence of an equivalent enrichment for paternally biased gene expression.Matrigenes that are favoured to be silenced in the context of the retinue behaviour would presumably be expressed from the patrigene.While we did not see an increase in paternally biased gene expression in our study, previous studies at a later point in worker bee behavioural maturation found evidence of higher paternally biased than maternally biased gene expression in both reproductive and brain tissues (Galbraith et al., 2016(Galbraith et al., , 2021)).Thus, there may be a temporal component to the silencing of matrigenes that is not being captured in the earlier timepoint of worker bee behavioural maturation under examination here.

| Genes showing intragenomic conflict are found within gene networks
Previous studies of the genetic basis of behavioural variation have revealed numerous insights into the coordination between genes within GRNs.Behavioural stimuli are often associated with predictable and reproducible changes in brain gene expression profiles (Alter et al., 2008;Benowitz et al., 2019;Friedman et al., 2020;Kabelik et al., 2021;Werkhoven et al., 2021;Whitfield et al., 2003;Wong et al., 2015), which are controlled by GRNs (MacNeil & Walhout, 2011).Thus, variation in the expression and regulation of genes within brain GRNs can be predictive of behavioural states (Sinha et al., 2020).In honey bee (Apis mellifera) colonies, female worker bees exhibit extraordinary examples of behavioural variation between nestmates (Smith et al., 2022) and are tractable models for genetic experiments (Toth & Zayed, 2021).Thus, studies in the honey bee have shown how transcriptional variation of genes within brain GRNs can contribute to variation in behaviour (Chandrasekaran et al., 2011;Ingram et al., 2011;Jones et al., 2020;Mikheyev & Linksvayer, 2015;Shpigler et al., 2017Shpigler et al., , 2019)).Here, we explored how conflicts within these genes-which are predicted to occur due to parent-specific differences in optimal strategies for maximizing reproductive success (Gardner & Úbeda, 2017)-can also contribute to behavioural variation.
Several studies have demonstrated that transcriptional variation of coregulated genes, versus large differences in individual genes, may be responsible for behavioural variation (Sinha et al., 2020).Additionally, meaningful transcriptomic differences associated with behavioural states may be missed by standard differential expression tests (Taylor et al., 2021).In our study, we did not detect any significantly differentially expressed genes between QMP-responsive and QMP-unresponsive bees, whereas a previous study of this behaviour identified hundreds of differentially expressed genes (Kocher et al., 2010).The behavioural states of the bees we selected for sequencing (from one reciprocal cross per Blocks 1, 3 and 4) were determined based on a single qualitative trait at a single timepoint rather than frequency of response to QMP across multiple timepoints (as in Blocks 5 and 6 in our study and in Kocher et al., 2010).
Thus, the transcriptomic differences associated with response to QMP in our study were likely too subtle to detect with standard differential expression tests.Therefore, we inferred co-expressed gene networks, identifying two modules that were correlated with behavioural state (Table 1; Figures S8 and S9).Additionally, we identified a subset of the transcriptome that was predictive of behavioural state in a SVM trained on the transcriptomic data (Figures S5 and   S7), and the genes in this SVM were enriched in a co-expression module (Table 1).
Our results provide support for the model of intragenomic conflict altering phenotypes by cascading impacts on downstream GRNs, versus directly altering expression of a single gene with major effects (Patten et al., 2016).As in previous studies in bees (Bresnahan et al., 2023;Galbraith et al., 2016Galbraith et al., , 2021;;Marshall et al., 2020), there was little overlap between genes that showed parent-biased expression and genes whose overall expression differences were associated with behavioural variation.However, our SVM and the combined set of genes showing parent-biased expression both showed overrepresentation for genes encoding proteins with DNA binding function, including several TFs.Moreover, we found that genes correlated with behavioural state, genes which predict response to QMP and genes engaged in conflict are enriched in conserved gene networks and functional pathways and are potentially regulated by a set of TFs engaged in intragenomic conflict.
One of the co-expression modules that was correlated with behavioural state (Module 7) contained genes with known functions in the retinue behaviour and behavioural maturation, including Vg and JHe.Oogenesis in insects involves the production of the egg yolk proteins from the egg yolk precursor Vg (Roy et al., 2018).
In honey bees, Vg and juvenile hormone (JH) are coregulated in a double-repressor relationship that regulates behavioural maturation in workers: High Vg titres suppress JH and foraging behaviour, while high JH levels suppress Vg and nursing behaviour (Harwood et al., 2017).In workers, QMP exposure suppresses JH biosynthesis, and JHe degrades JH (Mackert et al., 2008) delaying the transition from nursing to foraging (Pankiw et al., 1998).
The other co-expression module that was correlated with behavioural state (Module 4) was enriched for genes in the Notch signalling pathway.Chemical inhibition of Notch signalling in worker bees has been demonstrated to overcome the repressive effect of QMP on ovary activation (Duncan et al., 2016).Chemical inhibition of EZH2 (Tan et al., 2007), the enzymatic catalytic subunit of PRC2, has also been demonstrated to increase ovarian development in worker bees (Duncan et al., 2020).In our study, Module 4 was enriched for genes predicted to be regulated by the TF BtbVII-which was maternally biased in both QMP-responsive and unresponsive bees-including suppressor of hairless (Su(H)) and frizzled 2 (fz2).
Su(H) carries the Notch receptor signal to the nucleus where it acts as a transcriptional activator of enhancer of split complex genes (Falo-Sanjuan & Bray, 2020), suppressing neural development in D. melanogaster (Moore & Alexandre, 2020).BtbVII was demonstrated to physically interact with Krüppel homologue 1 (Kr-h1; Shokri et al., 2019), which is induced by JH and fz2 and regulates the expression of Vg (Jindra et al., 2015;Roy et al., 2018).Expression levels of Kr-h1 and fz2 in worker brains are downregulated by exposure to the queen or QMP (Grozinger & Robinson, 2007) and are higher in foragers compared with nurse bees (Fussnecker & Grozinger, 2008).
Interestingly, Module 7 was also enriched for genes predicted to be regulated by a Notch signalling inhibitor, lethal (2) giant larvae (Langevin et al., 2005), which was paternally biased in QMP-responsive workers.
One additional module was associated with response to QMP in that it was enriched for genes in our SVM classifier of behavioural state (Module 10).Both the SVM and module genes were overrepresented for structural constituents of chromatin and proteins with nucleosomal DNA binding function.This module was also enriched for genes predicted to be regulated by five parent-biased TFs, including fruitless (fru), scratch, stripe (referred to as early growth response protein 1 (Egr-1) in honey bees (Ugajin et al., 2013)), oestrogen-related receptor (ERR) and runt.Egr-1 has been implicated in the nurse-to-forager transition in worker bees, potentially by rapidly inducing neural development (Khamis et al., 2015;Lutz & Robinson, 2013;Shah et al., 2018).In D. melanogaster, both fru (Dalton et al., 2009) and ERR (Kovalenko et al., 2019) have been shown to regulate similar genes known to be induced by ecdysone receptor (EcR) signalling.In honey bees, EcR knockdown is associated with differences in transcript levels of several genes discussed above, including Kr-h1 (Mello et al., 2014), and QMP exposure is associated with high ecdysteroid titres (Trawinski & Fahrbach, 2018).In our study, fru was paternally biased in unresponsive workers, whereas ERR was maternally biased in responsive workers.

| Potential role for chromatin regulation in parent-of-origin intragenomic conflicts
Genes engaged in conflicts of origin 'disagree' on how to use maternal resources (Gardner & Úbeda, 2017).In mammals, this disagreement occurs within genes that regulate placental development and nutrient uptake by the fetus (Barlow & Bartolomei, 2014).In plants, this conflict has been observed in genes that regulate nutrient deposition and storage in the endosperm of seeds, flowering time, fruit development and root growth (Rodrigues & Zilberman, 2015).In bumble bees, parent-biased genes in ovaries are associated with female gamete generation, regulation of ovulation, response to pheromone and histone H3 acetylation (Marshall et al., 2020).Similarly, in honey bees, intragenomic conflict in worker ovary activation occurs in genes associated with neuron differentiation and embryonic development in fat body and ovaries (Galbraith et al., 2016) and reproduction-related hormonal signalling in brains (Galbraith et al., 2021).In our study, parent-biased transcripts in QMP-responsive bees were enriched for chromatin binding function, including several polycomb group genes (Pc, Sfmbt and RING1) and trithorax group genes (trx, Mi-2, polybromo, MTA3, Wdr82 and Sbf), suggesting a potential role for chromatin regulation in mediating intragenomic in the retinue behaviour.
Recent studies in plants and mammals have revealed a noncanonical epigenetic mechanism for genomic imprinting, which occurs through allelic differences in chromatin accessibility mediated by histone modifications, rather than allele-specific transcriptional repression by DNA methylation (DNAme; Andergassen et al., 2021;Batista & Köhler, 2020;Hanna & Kelsey, 2021).In insects, except for paternal genome elimination, genomic imprinting appears absent (Coolon et al., 2012;Pegoraro et al., 2017), and DNAme does not repress transcription (Duncan et al., 2022).Thus, in honey bees and bumblebees, allele-specific differences in DNAme are not associated with allele-specific differences in transcription (Marshall et al., 2020;Wu et al., 2020).However, a mechanism of allelespecific transcriptional regulation like noncanonical imprinting has been described in D. melanogaster (Floc'hlay et al., 2021).Given that polycomb repressive complexes (PRCs) are more widely conserved than DNAme among organisms (Schuettengruber et al., 2017), analogous mechanisms of noncanonical imprinting may be present

1
and 2, or C1 and C2) and two EHB colonies of A. m. ligustica stock (C3 and C4); these were managed at the Janice and John G. Thomas Honey Bee Facility of Texas A&M University (TAMU) in Bryan, TX.In studies conducted in 2021, EHB queens and drones were sourced from four colonies of A. m. ligustica 'Cordovan' stock (C5-C8; Koehnen & Sons, Inc., Glenn, CA), two colonies of a mixed A. m. ligustica stock (C9 and C10; BeeWeaver Honey Farm, Navasota, TX), one colony of A. m. carnica stock (C11) and one colony of A. m. caucasia stock (C12; Old Sol Apiaries, Rogue River, Oregon); these were managed at a Penn State University apiary in University Park, PA.The source colonies were separated into blocks, with two F0 colonies from different genetic backgrounds assigned to each block (Figure 1).
Three QMP-responsive and three QMP-unresponsive individuals of each maternal lineage (12 per 3 reciprocal crosses) were selected for sequencing.Heads were detached, and whole brains were dissected on dry ice and placed immediately in TRIzol (Invitrogen).Total RNA was isolated via TRIzol/chloroform extraction followed by precipitation with isopropanol.Samples were submitted to the Penn State Genomics Core for quantitation via Qubit and assessment of RNA integrity by RNA ScreenTape analysis on a TapeStation 4150 (Agilent, Santa Clara, CA).In total, 36 sequencing libraries were prepared using Illumina TruSeq Stranded mRNA Prep kits with unique 5′ and 3′ indexes, pooled and sequenced on a NextSeq 2000 P3 flow cell to generate approximately 34 million 150-bp single-end reads per sample (Table Parent-of-origin effects on the retinue response behaviour.(a) Blocks 1-4, the proportion of bees responding to QMP of each maternal lineage (source colony) from two reciprocal crosses per block.Significance of one-way ANOVAs is indicated (**p < .01).(b) Blocks 5 and 6, the mean frequency of response to QMP of individual bees, indicated by points, from six cage colonies constructed from each of three reciprocal crosses per block.Mixedeffects model significance of the effect of maternal lineage nested within blocks is indicated (**p < .01,***p < .001).

F
Figure 4 for results from Block 3 and Figures S3 and S4 for results from Blocks 1 and 4).Parent-biased transcripts in responsive bees were enriched for genes encoding proteins with chromatin binding function, including several core components of polycomb and trithorax group complexes.Overall, parent-biased genes were enriched for protein and sequence-specific DNA binding (TableS25).
Column 1: module name and size (number of genes).Column 2: the correlation coefficient and associated p-value for correlation with behavioural state.Green shading indicates significant correlations (p < .05).Columns 3-4: the number of genes in the module that overlapped with the given gene list and the associated p-value for these overlapping genes in Fisher's exact tests.Green shading indicates significant enrichment (p < .05).