The effect of maternal care on gene expression and DNA methylation in a subsocial bee

Abstract

Developmental plasticity describes the influence of environmental factors on phenotypic variation. An important mediator of developmental plasticity in many animals is parental care. Here, we examine the consequences of maternal care on offspring after the initial mass provisioning of brood in the small carpenter bee, Ceratina calcarata. Removal of the mother during larval development leads to increased aggression and avoidance in adulthood. This corresponds with changes in expression of over one thousand genes, alternative splicing of hundreds of genes, and significant changes to DNA methylation. We identify genes related to metabolic and neuronal functions that may influence developmental plasticity and aggression. We observe no genome-wide association between differential DNA methylation and differential gene expression or splicing, though indirect relationships may exist between these factors. Our results provide insight into the gene regulatory context of DNA methylation in insects and the molecular avenues through which variation in maternal care influences developmental plasticity.

Fig. 1

Behavioral response of adult offspring to early termination of maternal care. a Mothers establish nests in dead broken stems and produce all offspring in a linear hollow. Mothers typically remain on the nest to guard and groom immature offspring (care) or were removed from nests (no care). b Drawings of adult offspring engaging in aggression and avoidance behaviors in circle tube arenas. Illustrations in a and b courtesy of Wyatt Shell (used with permission). c Behaviors exhibited by offspring in circle tube assay in the presence and absence of maternal care (p-values are generated from two-tailed Mann–Whitney U-tests, n = 52; whiskers show observations within 1.5 × IQR of the lower and upper quartiles)

Fig. 2

Loss of maternal care leads to extensive changes in gene expression. a Hierarchical clustering of RPKM values (n = 11,068 genes) using the “complete” method with correlation-based distances indicates that global gene expression levels cluster into “No Care” (N) and “Care” (C) groups, with high support values using both the approximately unbiased (au) and bootstrap probability (bp) statistics. b Volcano plot showing the genes considered differentially expressed (blue) and those not significantly differentially expressed (red). Genes that had more than a twofold difference in gene expression are shown outside of the dashed lines. The numbers of differentially expressed genes are indicated above

Fig. 3

Loss of maternal care leads to localized, significant changes in DNA methylation. a Hierarchical clustering of gene level measures of DNA methylation (n = 19,361 genes) using the “complete” method with correlation-based distances reveals no global clustering by sample type (N, no care; C, care), with low support values using both the approximately unbiased (au) and bootstrap probability (bp) statistics. Thus, this tree can be considered one polytomy. Similar results are obtained when clustering only methylated genes. b Volcano plot showing the 200 bp regions considered differentially methylated (blue) and those not considered significant (red), as well as the thresholds used to define these categories (black lines). Methylation difference is given as a percentage. c Quadratic Assignment Protocol computed null distribution for the number of DMRs possibly occurring randomly using 50,000 randomized permutations and the associated number of experimental DMRs. The red line indicates the observed number of DMRs

Fig. 4

Overlap among differentially expressed, alternatively spliced, and differentially methylated genes. a UpSet plot illustrating the intersection between genes containing differentially methylated regions called by MethylKit, alternatively spliced genes called using DEXseq and rMATS, and differentially expressed genes computed using edgeR. The nature of a given intersection is indicated by the dots below the bar plot. For instance, the 99 genes in the fourth column are alternatively spliced and differentially expressed but not differentially methylated. b–d DEXseq plots showcasing the instances where a DMR directly overlapped a differentially utilized exon. The top panel shows average exon usage for each treatment group. The gene model is represented below each exon usage plot with purple exons indicating statistical support for differential usage (FDR-corrected p-value <0.01). Direction of transcription is indicated by an arrow and DMRs are overlaid beneath the gene model with colors indicating the change in methylation level