Changes of gene expression but not cytosine methylation are associated with male parental care reflecting behavioural state, social context and individual flexibility

Behaviour, like all phenotypes, is traceable to how and when genes are expressed. Transcriptional profiling has revealed distinct transcriptional architectures associated with distinct behavioural states (Zayed and Robinson, 2012; Cardoso et al., 2015; Parker et al., 2015; Palmer et al., 2016; Jacobs et al., 2016), which is further reflected in protein abundance (Cunningham et al., 2017). However, gene expression differences within and across behavioural states are not only a result of the genes that underpin the behaviour itself, but can also reflect other factors, such as social context and simple individual variation. An outstanding question is how much gene expression is associated with each source of variation when more than one factor is manipulated. Assessing the influence of multiple factors requires examining those factors with experimental designs that minimize other differences that can exist when examining highly distinct behavioural states (Benowitz et al., 2017; Kronauer and Libbrecht, 2018) and that allow a partitioning of effects. Furthermore, the mechanisms that regulate the variation of gene expression differences are not fully characterized, but need to be to fully understand of the evolution and mechanistic basis of behaviour.

In this study, we sought to estimate the association of gene expression with three different factors: (1) differences in behavioural state, (2) differences in social context and (3) differences in individual flexibility. We examined all three of these factors in terms of how they influence male parental care behaviour of the subsocial beetle Nicrophorus vespilloides Herbst 1783. This social behaviour displays considerable plasticity, making it accessible for the investigation of transcriptional architecture of flexible social behaviours under different conditions. In this species, males but not female parental care is naturally plastic (Smiseth et al., 2005). With a mate, males do not always (but can) participate in the feeding of the offspring, but instead provide indirect forms of care, such as excretion of anti-microbial compounds to cover the carcass (Smiseth et al., 2005; Capodeanu-Nägler et al., 2018). With the removal of his female mate, some males begin feeding offspring (Smiseth et al., 2005). Because we can manipulate male parental care via changes in social context, we can generate factorial crosses of males with or without mates that do or do not feed offspring. This helps us directly isolate associations of gene expression with behavioural state, social context and individual flexibility displayed for those behaviours (Table 1). This is important because social behaviour is multifaceted, and moving to factorial designs helps us to begin disentangling single variables, rather than comparing gene expression across behavioural states that necessarily differ by many variables (Lockett et al., 2012; Benowitz et al., 2017).

Within this experimental design, we also tested whether cytosine methylation regulated any rapid changes in gene expression during socially responsive parental care. Cytosine methylation is a core mechanism regulating gene expression (Cardoso et al., 2015). It is stable (Turecki, 2014; Yan et al., 2015), reversible (Yan et al., 2015) and can have relatively short-term turnover in animals (Levenson et al., 2006; Guo et al., 2011; Herb et al., 2012; Mizuno et al., 2012; Baker-Andresen et al., 2013; a contrario sensu Cardoso et al., 2015). Therefore, cytosine methylation has excellent characteristics to regulate the gene expression of underlying behaviours (Cardoso et al., 2015; Yan et al., 2015). Cytosine methylation is associated with different behaviours of a range of insects, including different hymenopterans (Kucharski et al., 2008; Lyko et al., 2010; Bonasio et al., 2012; Foret et al., 2012; Herb et al., 2012; Lockett et al., 2012; Amarasinghe et al., 2014; Kucharski et al., 2016; Herb et al., 2018) and an orthopteran (Wang et al., 2014). However, its association with behaviour is not ubiquitous, as cytosine methylation is not associated with different behaviours of several hymenopterans (Patalano et al., 2015; Libbrecht et al., 2016; Glastad et al., 2017; Toth and Rehan, 2017), nor with the evolution of social behaviour of insects in general (Bewick et al., 2017). We have a sequenced and annotated genome for N. vespilloides, and there is gene body cytosine methylation allowing us to directly assess the association of cytosine methylation with different factors (Cunningham et al., 2015). The role of cytosine methylation underlying gene expression differences of transient behaviours has not been assessed. More generally, it is still unknown whether the mechanisms underlying rapid, transient and flexible transitions of behaviour are the same as those that are associated with longer-term behavioural transitions.

Our first goal was to identify gene expression associated with three different factors: differences in behavioural state (directly caring versus non-caring), differences in social context (with or without a mate) and differences in individual flexibility (whether an individual changed their behaviour during the study; Table 1). We predicted that differences in behaviour would have a large influence on gene expression (Parker et al., 2015), followed by differences in social context (Parker et al., 2015), and the influence of individual flexibility of behaviour is largely unknown. We also predicted the possible change in expression of several pathways, including neuropeptides (Cunningham et al., 2016, 2017; Bukhari et al., 2017), neural remodelling factors (e.g. bdnf; Cunha et al., 2010) and genes associated with transcriptional regulation in general (Cardoso et al., 2015). Our second goal was to assess whether cytosine methylation underpinned the rapid changes in gene expression seen during rapid changes in behaviour using whole genome bisulfite sequencing (WGBS) of the same males used for the gene expression experiments. Assuming cytosine methylation underpins behavioural changes, we expected that cytosine methylation levels would change for behaviourally responsive genes (Cardoso et al., 2015; Yan et al., 2015). We found many differences in gene expression between caring and non-caring behavioural states, and fewer expression differences owing to social context or individual flexibility. Very few cytosine methylation changes were associated with any of the factors we manipulated. Thus, differential expression of genes accompanies rapidly changing behaviour with a hierarchy from behavioural state, social context and individual flexibility; however, cytosine methylation does not appear to underpin any of these rapid changes and the epigenetic mechanisms that influence this process remain to be identified.

The parental care behaviour of N. vespilloides is multifaceted, easily observed and reliably scored (Smiseth and Moore, 2004; Smiseth et al., 2005; Walling et al., 2008). Parental care in all burying beetle species is extensive and elaborate, including direct provisioning of regurgitated food to begging offspring (Eggert and Müller, 1997; Scott, 1998). Parental care can be uniparental or biparental, often within a species. Adults search for and bury a small vertebrate carcass on which they feed and rear offspring. Parents provide both indirect and direct care. Before offspring are born there is indirect care involving stripping the fur (or feathers or scales) from the carcass, forming it into a nest, and preventing microbial growth on the carcass through excretions (e.g. Palmer et al., 2016; Duarte et al., 2018; Shukla et al., 2018; Wang and Rozen, 2018). The latter form of indirect care also occurs after offspring are present, along with resource defence (Walling et al., 2008). Eggs are deposited away from the carcass while it is being manipulated into a suitable larval food resource. When eggs hatch, the larvae crawl to the carcass and reside in a small cavity excised by the parents in the carcass. Parents provide direct care by regurgitating pre-digested carrion directly to their dependent, begging offspring and by depositing oral secretions into the larval cavity to prepare it for larval occupation and consumption (Scott, 1998). In N. vespilloides, the species studied here, parental care can be provided under multiple social contexts; by both parents or either individually without influencing the survival or vigour of larvae (Parker et al., 2015). When both parents are present, females bias their time toward direct care of offspring (i.e. feeding them) while males spend more time on indirect care (Smiseth et al., 2005). This system is amenable to our experimental manipulation (Table 1), as removing females while larvae are still young results in some males initiating direct care (Smiseth et al., 2005), while other males will maintain their level of direct care observed when a female is present (generally, high level or none). This allowed us to create a factorial design to separate the influence of behavioural state, social context and individual flexibility on gene expression.

We obtained beetles from an outbred colony of N. vespilloides, originating from Cornwall, UK, and maintained at the University of Georgia (Cunningham et al., 2014, 2017). This colony is augmented with new families yearly from the same origin population. We followed the protocol of Smiseth et al. (2005) to generate flexibly caring males. Unrelated female and male pairs (age 14–29 days) were placed into a mating box with a mouse carcass (19–21 g) and 2.54 cm of moist soil. The boxes were observed every morning (approximately 09:00 h) and evening (approximately 17:00 h) starting at 60 h post-pairing until larvae arrived at the carcass. Twenty-one hours after larval arrival, each pair was observed using 1 min scans for 10 min per hour for four observation periods. We then repeated the observation protocol 24 h later. There were two treatments on day 2: we removed the female from pairs where the males showed no direct care on day 1, and left the pairs intact for the other half. If males were observed caring on day 1, we left the pair intact. All pairs were observed both days regardless of treatment.

Because we were first interested in separating the influence of three factors on gene expression, we designed an experiment that manipulated males into one of four different treatments that allowed us to assess three a priori contrasts (Table 1). The first sample, phenotypically ‘flexible care’, contained males that initiated direct care (directly feeding offspring) when the female was removed. The second sample, ‘non-flexible no-care’, contained males that were not observed to directly care for offspring even if the female was removed. The third sample, ‘biparental no-care’, contained males that never directly care for offspring with the female present both days. The fourth sample, ‘biparental care’, contained males that that always directly cared for offspring with the female present both days. No-care samples might have dispersed from the carcass, which is why they might not have been observed to care. It is hard to tell with our protocol, as the boxes are so small that males are often found near but not on the carcass. This can signify rest or dispersal, which can only be distinguished with multiple widely spaced observations that directly opposes our ability to gather samples directly after the behavioural observations. Other samples are not available as males do not provide direct care on day 2 if they do not provide direct care on day 1 in the presence of females, and males that provide direct care on day 1 rarely change to no direct care on day 2 regardless of the presence or absence of the female. To maximize power, we only selected males for analysis that showed ‘pure’ phenotypes; that is, consistently high care or absolutely no direct care throughout all observation periods. We observed 104 pairs in total (Table S1).

We dissected brains from individual males as in Cunningham et al. (2014), with the exception that samples were snap-frozen in liquid nitrogen after dissection and stored at −80°C. From these samples, we extracted RNA and genomic DNA (gDNA) simultaneously using Qiagen's AllPrep DNA/RNA Mini Kit (cat. no. 80284; Hilden, Germany) following the manufacturer's protocol after homogenization with a Kontes handheld pestle (Kimble Chase, Rockwood, TN, USA) to allow us to quantify both gene expression and methylation level from the same individual. We quantified RNA and gDNA with a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) using the RNA Broad Range and dsDNA High Sensitivity protocols, respectively, following the manufacturer's instructions.

We prepared libraries for RNA-seq with a modified Smart-seq2 protocol (Picelli et al., 2014) using a target of 80 ng of total RNA per library and barcoded with Illumina TruSeq indices. Libraries were subjected to solid phase reversible immobilization (SPRI) to select for fragments between 300 and 1000 bp, and insert size was estimated with a Fragment Analyzer Automated CE System (Advanced Analytical, Ankeny, IA, USA). We sequenced 24 samples (six from each of the four behavioural states), assigned to one of two pools to evenly distribute samples based on experimental factors across the two lanes, with a 75 bp single-end (SE) protocol using Illumina's NextSeq500 with a high-output flow cell targeting 35 million reads per sample at the University of Georgia's Georgia Genomics Facility (Table S2).

We assessed the quality of the raw sequencing reads using FastQC (v0.11.4; default settings; https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We trimmed for the transposase adapter, with reads based on quality (Phred>15 at both ends), trimmed the last two base pairs of the reads owing to highly skewed nucleotide frequencies, and reassessed the quality of the reads using FastQC (v0.11.4) and Cutadapt (v1.9.dev1; Martin, 2011).

We combined data from the four groups of males with different experiences of parental care, social context and individual flexibility to parse the effect of different influences and refine the potential causal differences (Table 1). We performed three contrasts. First, we compared all those individuals that displayed direct care with those that did not, regardless of social context (behavioural state contrast). This compares males that transition from a no-care state to a direct care state with those that do not make this transition with or without the female. Second, we compared individuals in the presence of a female both days with those where a female was absent the second day, regardless of their own behaviour (social context contrast). This tests for the transcriptional response to a difference in social context. Third, we tested for the transcriptional response to changes in individual behaviour regardless of social context or starting behaviour (individual flexibility contrast). This directly compares flexible with non-flexible individuals. Taken together, these three contrasts lead to a description of gene expression unique to each source of variation influencing behavioural changes.

We used RSEM (v1.2.20; default settings; Li and Dewey, 2011) with BowTie2 (v2.2.9; default settings; Langmead and Salzberg, 2012) to map and quantify reads against the N. vespilloides official gene set Nv OGS v1.0 transcriptome (Cunningham et al., 2015). To better assess the completeness of the Nv OGS v1.0 before mapping, we used the updated BUSCO gene set (v2.0; default settings; Simão et al., 2015) with the Arthropoda hidden Markov models (2675 HMM gene models). This gene set is defined as gene models that are present in 90% of the searched species as single-copy orthologues. We found that 2607 (97.5%) genes were present with 2484 (92.8%) as complete genes. Of the complete genes, 2183 were single-copy orthologues and 301 were duplicated. A further 123 genes were fragmented.

Differential expression was estimated using both parametric and non-parametric differential gene expression analysis to find genes that individually exhibited strong responses to our manipulation. These two methods find differential expression based on different biological signals and so can identify different sets of genes between contrasts of interest. For each analysis, we performed three contrasts (Table 1).

We imported the expected read count per gene from RSEM into the DESeq2 package (v1.12.3; default settings; Love et al., 2014) using the tximport package (v1.0.3; Soneson et al., 2015) of R (v3.3.1; https://www.r-project.org/). Following the suggested workflow of DESeq2, we performed overall sample quality control by visual inspection of a principal component analysis (PCA) plot of the data and removed two outlier samples (one flexible care, one non-flexible care). Statistical significance was assessed after a Benjamini–Hochberg (BH) correction of P-values (Benjamini and Hochberg, 1995). We used NOISeq (2.16.0; Tarazona et al., 2015) to test for differential gene expression as a non-parametric method. Following the suggested workflow of NOISeq, we performed overall sample quality control by visual inspection of a PCA plot using trimmed mean of M-values (TMM) standardized data and removed one outlier sample (flexible care), as per program guidelines. Each analysis was conducted using TMM standardized data, filtering out genes with counts per million reads (CPM) <1, correcting for gene length, substituting zero gene counts with 0.01 to avoid undefined gene counts, and with 20 permutations using the NOISeqBIO function. Statistical significance was assessed after a BH correction of P-values. We used all the unique differentially expressed genes from either the DESeq2 or NOISeq gene list (i.e. the union of the two gene lists) for each of the three contrasts to test for enrichments of all three categories of gene ontology (GO) terms: biological process, molecular function and cellular component. We used the AgriGO webserver to test for enriched GO terms (Du et al., 2010). We performed a singular enrichment analysis (SEA) using complete GO terms and a hypergeometric test with a BH correction. The complete list of GO terms assigned to all N. vespilloides genes was used as the background for the enrichment test.

Because genes usually act within a network, and whole networks can exhibit responses to a manipulation even if the individual genes within the network do not, we also performed a weighted gene co-expression network analysis (WGCNA). This technique also allows for the centrality of a gene to a network to be estimated with the assumption that genes deeply connected within a network are of increased overall importance because changing their expression influences many other genes. We again looked for associations with our three contrasts and the expression of gene modules between these contrasts. We used the WGCNA package of R (Storey, 2002; Langfelder and Horvath, 2008) to perform a WGCNA using default guidelines and parameters. We used the variance stabilized transformation, which is blind to the study design from DESeq2, with the same two outlier samples removed, as input data, with genes with <10 reads in 20 samples removed as per the program’s suggestion. We converted the correlation matrix of variance stabilized transformed values (DESeq2's default transformation) to a signed adjacency matrix with an exponent of 10 and a minimum module size of 30. We tested for an association between modules and traits of interest using the biweight mid-correlation (bicor) function with a robustY setting, as per program guidelines for our data types. Modules significantly associated with traits were assessed for enrichment of GO terms as described above.

We also established orthology between Amel OGS3.2 (Elsik et al., 2014) and Nicve OGS1.0 using OrthoFinder (default setting using diamond; v2.1.2; Emms and Kelly, 2015) to assess the overlap of differentially expressed ‘direct caring’ genes between Shpigler et al. (2018) and the present study.

We used MethylC-seq to estimate levels of cytosine methylation associated with different behavioural states. We prepared MethylC-seq libraries following Urich et al. (2015) targeting 200 ng of gDNA as input per library. Six individuals, three each from sample groups 1 and 2 (Table 1), that we used for RNA-seq were haphazardly chosen for whole genome bisulfite sequencing. Libraries were quality controlled with the above RNA-seq protocol. We sequenced six adult samples with a 150 bp SE protocol using Illumina's NextSeq500 with a high-output flow cell at the University of Georgia's Georgia Genomics Facility (Table S2).

We followed the protocol of Cunningham et al. (2015) to determine the methylation status of individual cytosines and genes that were used to survey the methylome of larval N. vespilloides. Briefly, we used the methylpy analysis pipeline (Schultz et al., 2015), which checks reads for adapter contamination and quality score trimming with cutadapt (v1.9dev), maps with Bowtie1 (v1.1.1; parameters: -S -k 1 -m 1 --chunkmbs 3072 --best --strata -o 4 -e 80 -l 20 -n 0), removes PCR duplicate reads with Picard (v2.4.1; default settings; broadinstitute.github.io/picard), and uses a BH-corrected binomial test against the sample specific non-conversion rate of fully unmethylated lambda gDNA to call methylated cytosines. Cytosines within a region of interest [here, coding sequence (CDS)] were aggregated and a BH-corrected binomial test against the mean percentage of methylated cytosines per gene is used to call methylated genes. To estimate how conserved gene methylation status is between adult and larval life history stages, we re-analyzed the six adult samples from this study and the three larvae samples from Cunningham et al. (2015; NCBI BioProject: PRJNA283826) together. To address the influence of different sequencing coverage between these samples, we restricted our analysis to genes that had at least five CpGs covered with at least three mapped reads within the CDS regions for all nine samples (i.e. we only assessed genes with sufficient amounts of information from all samples to reduce the influence of noise from low-coverage CpGs and coverage differences between samples) (Cunningham et al., 2015). A BH-corrected binomial test determined the methylation status of each gene within each sample using the mean percent of methylated CpGs of all samples across all genes as the null probability. Genes identified as methylated in all adult samples and unmethylated in all larval samples were defined as adult-specific methylated genes, and vice versa. We defined the overlap as the union of adult methylated genes compared with the union of the larval methylated genes.

We estimated differential cytosine methylation amongst the two adult behavioural states (flexible care versus nonflexible no-care) in two different ways (qualitative and quantitative) at the gene (Patalano et al., 2015) and individual nucleotide (Libbrecht et al., 2016) levels. Our analysis was designed within an exploratory framework to capture any signal of individual cytosine or gene methylation status associated with social behaviour. For the qualitative analysis at the gene level, we assessed how many genes were consistently methylated or non-methylated in one sample group while having the opposite methylation status in the other sample group. The quantitative analysis was a BH-corrected t-test of the proportion of methylated cytosines across a gene or a BH-corrected t-test of weighted methylation level across a gene (number of methylated reads/all reads mapped to a cytosine; Schultz et al., 2012) with at least 10 mapped cytosines (12,627 genes meet the minimum coverage threshold; Patalano et al., 2015).

For the qualitative analysis at the nucleotide level, we assessed how many cytosines were methylated or non-methylated in one sample group while having the opposite methylation status in the other sample group. The quantitative analysis was a BH-corrected t-test of the weighted methylation level (number of methylated reads/all reads mapped to a cytosine) for every cytosine that was mapped in all adult samples with at least five reads.

In the sample in which males were induced to shift from no care to direct care (Table 1, sample 1), the percentage of observed time spent directly feeding larvae shifted from 0% (with female; day 1) to 28.3±0.4% (after female removal; day 2). In samples where females were not removed but male direct care was observed (Table 1, sample 4), males spent 34.0±5.5% of the observation period on direct care in day 1, and 35.9±4.1% of the observation period directly caring for larvae on day 2. These results corroborate those of Smiseth et al. (2005).

There was no statistically significant difference in the proportion of direct care (DC) and indirect care (IC) for males that expressed care with and without a partner (day 1: DC–34.0%, IC–4.6%; day 2: DC–35.9%, IC–3.5; DC: t_23.5=−0.22, P=0.826; IC: t₂₄=0.269, P=0.79). There was also no statistically significant difference in the proportion of DC and IC for males that expressed care (day 2: biparental care–DC–35.9%, IC–0.35%; flexible care–DC–28.3%, IC–5.3%; DC: t_25.2=1.010, P=0.32; IC: t_24.2=−0.522, P=0.61).

To identify differentially expressed genes and gene co-expression networks associated with changes of behaviour in different contexts, we investigated gene expression between three contrasts: behavioural state, social context and individual flexibility (Table 1; Table S3). For the behavioural state contrast, we found 522 differentially expressed genes using parametric analysis (177 genes upregulated and 345 genes downregulated contrasting caring to non-caring), 150 differentially expressed genes using non-parametric analysis (138 genes upregulated and 12 genes downregulated; the union of the DESeq2 and NOISeq gene sets is 552 genes; Fig. 1) and seven co-expressed gene modules using WGCNA (modules 1, 2, 5, 7, 8, 9 and 10; Table 2). For the social context contrast, we found 97 differentially expressed genes using parametric analysis (45 genes upregulated and 52 genes downregulated comparing without and with partner), zero genes differentially expressed using non-parametric analysis and one co-expressed gene module using WGCNA (module 5; Table 2). For the individual flexibility contrast, we found 17 differentially expressed genes using parametric analysis (eight genes upregulated and nine genes downregulated comparing flexible with non-flexible), three differentially expressed genes using non-parametric analysis (all upregulated; union of two sets is 19 genes) and three co-expressed gene modules using WGCNA (modules 1, 9 and 10; Table 2). A high number of downregulated genes is consistent with the pattern of gene expression sampled during caring and compared with the post-care period (Parker et al., 2015). As expected, there was little overlap between the differentially expressed genes between the contrasts, suggesting that we could cleanly dissect each effect (Fig. 2).

We next used GO analysis to examine the potential functions or functional categories of the genes and gene co-expression networks associated with each contrast. We found 77 GO terms enriched for the behavioural state contrast, with glutamine family amino acid metabolism, cellular aromatic compound metabolism, carboxylic acid metabolism, oxoacid metabolism, cellular amino acid biosynthetic processes and organic acid metabolism being the most significantly associated (all P=0.0063; Table S4). Only two of the seven gene co-expression networks associated with the behavioural state contrast had significant GO enrichment. Module 7 was enriched for terms related to mitochondria, cell envelope and organelle envelope (all P=0.037), whereas module 9 was enriched for terms related to cellular amino acid metabolism, carboxylic acid metabolism, oxoacid metabolism, organic acid metabolism and small molecule biosynthetic processes (all P=0.019). Genes differentially expressed associated with variation owing to difference in social context were enriched for GO terms related to only three terms: ion binding, cation binding and metal-ion binding (all P=0.011). The one gene co-expression network associated with social context had no significant GO enrichments. The differentially expressed genes of the individual flexibility contrast were not enriched for any GO terms. Of the three gene co-expression networks associated with the individual flexibility contrast, only module 9 had enriched GO terms (see above).

We investigated differences in gene or cytosine methylation to assess their relationship with flexibility in expressing care, focusing on a comparison of individuals that changed from no care to care and those that never changed (Table 1; sample 1 versus sample 2). This comparison should capture any mechanism associated with changes in behavioural state or individual flexibility. The genes methylated in reproductive adults overlapped highly with methylated genes in N. vespilloides larvae (99.4%; Fig. 3). However, we found that only 2.1% of conserved adult methylated genes were also differentially expressed in any of our three contrasts (Fig. 4, showing largest overlap contrast; Table S5).

We next asked whether any methylation changes at the gene level were associated with individual flexibility of adults. We found no association between the total number of methylated genes and changes in behaviour (t₄=0.714, P=0.515). We then asked whether methylation of individual genes differentiates these samples. We found 17 genes displaying a qualitative difference in methylation status. However, two methods of quantitative gene methylation analysis – percent of methylated cytosines and weighted methylation level – showed that zero and one gene, respectively, differed between flexibly expressed care and non-flexible no-care males.

It could be possible that methylation differences of individual cytosines (rather than across the entire gene body) are responsible for producing phenotypic differences. Therefore, we examined whether methylation of individual cytosines was associated with flexibly expressed care. Qualitatively, we found 460 cytosines with differing methylation status between the two groups. A permutation analysis of our samples showed that 510.5±307.0 (mean±s.d.) cytosines differed in methylation status. Therefore, 460 cytosines are no more than expected by chance, and provide little evidence that individual cytosine methylation is associated with behavioural state or individual flexibility. Furthermore, quantitative analysis of cytosine-weighted methylation level showed only a single nucleotide (out of 56,753 methylated cytosines that had coverage in all samples) significantly associated with behavioural state or individual flexibility.

Our results suggest a hierarchy of associations with gene expression during socially responsive parental care. Greater differences of gene expression were induced by manipulating behavioural states (directly caring versus non-caring), fewer associated with social context, and the least associated with individual variation in expressing a behaviour. The first result is consistent with the large body of studies showing differences between many behavioural states are strongly associated with gene expression differences and to a lesser extent with other factors (Zayed and Robinson, 2012; Cardoso et al., 2015; Parker et al., 2015; Toth and Rehan, 2017; Tripp et al., 2018). Better understanding of why some factors are more impactful to transcriptional architecture will require more studies using factorial designs to discover general patterns. The more flexible, and therefore rapid, the behavioural change, the fewer gene expression changes involved. This suggests that increased flexibility requires fewer changes of gene expression, perhaps because fewer phenotypic components are changing.

By going beyond broad state comparisons, we show that transcriptional architecture depends on the context in which behaviour is expressed. Most studies that compare differences in caring and non-caring states have compared different castes or developmental states in Hymenoptera (Toth and Rehan, 2017). Such transitions are typically long-lasting, and can include confounds from developmental change. It is unclear how to compare such studies and the present study, which is more concerned with rapid changes and switches that can be reversed. A more comparable study is that of Shpigler et al. (2018), which compared young nurses present with queen larvae to initiate alloparental care. Shpigler et al. (2018) sampled the mushroom body of the brain of nurses at 30, 60 and 120 min after care was initiated followed by RNA-seq. We examined the overlap of the genes identified by Shpigler et al. (2018) by establishing one-to-one orthology and then investigated the overlap between our 552 DEGs for the caring versus non-caring contrast and their 1268 DEGs from any sampling point. We found 49 orthology groups that contained genes that were differentially expressed in both gene sets (Table S6). We suggest this small overlap has to do with the difference of the sampling points, although it may also reflect the restricted brain region sampled in Shpigler et al. (2018). The latter hypothesis seems less likely, however, because our samples included mushroom bodies as well as other parts of the brain. Instead, we may have missed the genes causing a ‘change of behaviour’ and sampled the ‘stable expression of behaviour’ genes. We are clearly at a point where we need more functional studies of genes associated with parenting. This further highlights the message that considering the context, rapidity and type of change is important.

Comparisons can also be made with recent studies using N. vespilloides. Genes previous found to be differentially expressed in N. vespilloides when transitioning from a completely non-parenting to parenting state – such as vitellogenin 1, vitellogenin 2 and their receptor (vg1, vg2 and vgr; Roy-Zokan et al., 2015), and neuropeptide f and its receptor, npfr (Cunningham et al., 2016) – were not involved in the subtler plasticity examined here. Benowitz et al. (2017) sampled N. vespilloides and Nicrophorus orbicollis males and females that contrasted the transcriptional architecture of high versus low direct care expressed by parents in a 2×2×2 design. Their study found very subtle differences of transcription associate with the level of direct care expressed by parents. In the present study, we found very clear differences of gene expression when comparing males providing direct care with those that did not. In this way, we align much more with an earlier study of N. vespilloides by Parker et al. (2015) that compared very large differences in behaviours/social context (virgin to direct care to post-care) and found very clear differences in gene expression. The results of the present study reaffirm the suggestion of Benowitz et al. (2017) that transcription variation underpinning the phenotypic variation within a behaviour (high versus low direct care) is much subtler than the transcriptional variation that underpins complete changes of behaviours (high versus no direct care). We also suggest that this highlights how impactful small changes in behaviour can be to transcriptional architecture.

When we assessed the functional categorization of the differentially expressed genes and gene co-expression networks, we found an abundance of metabolic-related categories. Despite the abundance of GO terms related to metabolism, we do not expect these genes to reflect the energetic cost of parenting because we only sampled brains. Instead, we suggest that metabolic genes might be co-opted for a social function in N. vespilloides, as is argued elsewhere (Zayed and Robinson, 2012; Rittschof et al., 2014; Wu et al., 2014; Cunningham et al., 2016; Fischer and O'Connell, 2017). Alternatively, metabolic genes may be involved in neurotransmitter synthesis (Livingstone and Tempel, 1983), as many neurotransmitter pathways influence parental care (Mileva-Seitz et al., 2016). One potentially interesting candidate gene found in both the list of differentially expressed genes and as a hub gene in the gene co-expression network (module 9) associated with caring is NK homeobox 7 (nk7). This gene was also one of the only genes showing evidence for positive selection in the N. vespilloides genome (Cunningham et al., 2015), and thus multiple lines of evidence suggest it is an important regulator of parental care behaviour. The differentially expressed genes associated with differences in social context are related to ion binding, which might be associated with ion-gated channels in the brain that modulate neural activity (Simms and Zamponi, 2014). Thus, these channels may represent a candidate pathway mediating effects of the social context on behaviour. Individual flexibility of behaviour produced a clear gene expression signal, but the types of gene underlying this phenotype are difficult to classify. The gene co-expression network associated with flexibility is more strongly associated with caring than with individual flexibility per se. Individual flexibility in ants and bees is associated with morphological changes in the brain (Gronenberg et al., 1995; Groh et al., 2006), and thus we expected to detect genes annotated with neurotropic activity or neuron axon manipulation. That we made no such observation suggests that gross morphological changes in the brain might only be seen in species that make permanent or developmental changes between behavioural states (Cardoso et al., 2015). It is also possible that we sampled males too late to capture the genes involved in causal changes in behavioural change, especially the immediate early genes that respond within minutes to hours to a stimulus (Cardoso et al., 2015).

There is little evidence to suggest that methylation at the individual gene or individual cytosine level is associated with behavioural state or individual flexibility of male direct care in N. vespilloides. There were few differences at the gene or individual cytosine level between the two samples compared (flexible caring versus non-flexible non-caring). Furthermore, very few (2.1%) of the adult methylated genes were also genes that were differentially expressed for any of the three contrasts of gene expression. Adult methylated genes were also highly overlapping with larval methylated genes, which indicates that gene methylation is stable across broad life history stages (and generations), encompassing widespread behavioural and physiological changes.

Our results are informative because we assessed a transient behaviour, extending the range of behaviours that cytosine methylation has been thought to influence. Our results also fall in line with those of other studies of social insects demonstrating few differences of cytosine methylation between different behavioural states (Patalano et al., 2015; Libbrecht et al., 2016). Not all social Hymenoptera even have active DNA methylation systems (Standage et al., 2016). Therefore, cytosine methylation does not appear to be a general mechanism to regulate behavioural changes of insects (Patalano et al., 2015; Libbrecht et al., 2016; Bewick et al., 2017; Glastad et al., 2017), but it remains possible that it might regulate socially responsive gene expression of any one species, such as long-term changes of aggression after territorial challenges in honey bees (Herb et al., 2018).

Using the socially responsive and naturally variable male parental care of the subsocial beetle N. vespilloides, we made a series of comparisons to understand the influence of behavioural states, social context and individual flexibility on transcriptional architecture of a transient social behaviour. We found clear signals of gene expression after manipulating behavioural state (directly caring versus non-caring), social context (with or without a female mate) and, to a much lesser extent, an individual's ability to rapidly change behaviour. This suggests a complex and hierarchical influence on the transcriptional architecture of parenting behaviour by males. Research on behavioural transitions has long examined the role of single molecules, such as neuropeptides and hormones. Thus, it is perhaps unsurprising that an individual's ability to change behaviour might involve fewer changes of gene expression. Although differential gene expression has, for many genes, long been associated with changes of long-term or permanent behaviour (Zayed and Robinson, 2012; Cardoso et al., 2015), the present study further supports that gene expression is also associated with more rapid changes of behaviour responding to immediate environments. Contrary to previous predictions regarding epigenetic control of behaviour (Cardoso et al., 2015; Yan et al., 2015), we found no support for an association between cytosine methylation and the expression of direct parental care or individual flexibility, and conclude that rapid changes of cytosine methylation is not the mechanism underpinning the rapid behavioural transitions for this species. This leads to the conclusion that, contrary to some predictions, rapid gene expression affecting behaviour may be regulated by standard processes of transcriptional control. Our work suggests that studies of the genetic underpinnings of behavioural flexibility, a key attribute that defines behaviour as a unique phenotype (Bailey et al., 2018), should consider the extent to which the behavioural change is transient.

This is the 11th paper of the burying beetle Nicrophorus vespilloides molecular genetics of parenting project at the University of Georgia. We would like to thank Dr Sandra Steiger for thoughtful comments that improved this manuscript. We thank the University of Georgia's Georgia Advanced Computing Resource Center for computational infrastructure and technical support.