Parental care shapes the evolution of molecular genetic variation

Abstract Cooperative social behaviors, such as parental care, have long been hypothesized to relax selection leading to the accumulation of genetic variation in populations. Although the idea has been discussed for decades, there has been relatively little experimental work to investigate how social behavior contributes to genetic variation in populations. Here, we investigate how parental care can shape molecular genetic variation in the subsocial insect, Nicrophorus vespilloides. Using whole-genome sequencing of populations that had evolved in the presence or absence of parental care for 30 generations, we show that parental care maintains levels of standing genetic variation. In contrast, under a harsh environment without care, strong directional selection caused a reduction in genetic variation. Furthermore, we show that adaptation to the loss of care is associated with genetic divergence between populations at loci related to stress, morphological development, and transcriptional regulation. These data reveal how social behavior is linked to the genetic processes that shape and maintain genetic diversity within populations, and provides rare empirical evidence for an old hypothesis.


Introduction
While much recent work has focused on identifying genes that drive social behaviors (Cunningham, 2020;Toth & Robinson, 2007), relatively few studies have examined the longstanding hypothesis that social behavior affects the accumulation and maintenance of genetic variation.Yet, social living is associated with largescale restructuring and the evolution of genome organization and architecture (Rubenstein et al., 2019).In humans, benevolent social activities, such as modern health care, are thought to have led to the accumulation of deleterious mutations within populations (Kondrashov, 2017;Lynch, 2016).Therefore, the extent to which genetic variation is shaped by social behavior has implications for the health of populations and their capacity to rapidly adapt to environmental perturbations.However, there have been few empirical tests of how social behavior might drive genetic variation in practice.Here we investigate how a cooperative social behavior, namely the supply of parental care, contributes to genome-wide levels of genetic variation.We focus on parental care in a subsocial pair-breeding insect, rather than more elaborate forms of sociality, to avoid the confounding effects of extreme reproductive skew on genetic variation, which is common in cooperative insect societies (Social Insects, 1979).
Cooperative social interactions often function to shield social partners from a harsh physical environment, and the same is true for parental care (Cornwallis et al., 2017;Royle et al., 2012).Without cooperation generally, and care specifically, individuals would be exposed to strong, frequently directional, selection pressures from the abiotic environment, which would favor the evolution of new adaptations and cause an associated reduction in genetic variation.On the other hand, the presence of parental care relaxes selection from this wider environment, theoretically allowing genetic variation to accumulate.Indeed, several lines of evidence suggest that cooperative social behaviors, including care, can relax selection sufficiently to allow mildly deleterious mutations to accumulate within populations (Linksvayer & Wade, 2009;Pascoal et al., 2023;Pilakouta et al., 2015;Schrader et al., 2018;Snell-Rood et al., 2016).In this way, parental care could shift the "mutation-selection" balance by relaxing selection and preventing the elimination of new spontaneous mutations.The resulting increase in genetic variation could emerge in the form of single-nucleotide polymorphisms (SNPs) and/or other structural genetic variants (e.g., indels, transposable elements, and/or insertions/deletions) depending on the natural mutation rate of such variants.Exactly how care might maintain such variants has been the subject of some speculation (Snell-Rood et al., 2016).One possibility is that "cryptic" variants could be maintained in the population with a combination of care-induced genetic capacitors, epigenetic modifications, and/or RNA-mediated signals (Paaby & Rockman, 2014).Nevertheless, although the suggestion that cooperative social behavior can shape genetic variation is relatively longstanding, we still have a poor understanding of how and where it might cause change at a molecular genetic level.
Here, we use evolving populations of burying beetles (Nicrophorus vespilloides) to explore how parental care affects levels of standing genetic variation and how populations may adapt in the face of its loss.In natural populations of this locally abundant subsocial insect, burying beetle parents raise their young on a carrion nest, formed from a small dead animal, such as a mouse or songbird.There is continuous variation in the level of parental care supplied, with around 5% of parents abandoning the brood before their young have even hatched (Scott, 1998).Offspring can survive without parental care, at least in the laboratory.
We exploited this natural variation in care to establish two types of experimentally evolving populations in the laboratory, which varied only in the family environment that larvae experienced during development and where the same family environment was created for successive generations within populations.In Full Care populations (FC), parents remained with their young throughout development, whereas in No Care populations (NC), parents were removed just prior to hatching.No Care populations rapidly adapted to this regime (within 14 generations), with adaptive change being detectable through increases in breeding success and larval density (see Schrader et al., 2017).Moreover, we have previously shown that No Care populations evolved adaptively (Schrader et al., 2015) and divergently from Full Care populations in the extent of the prehatching care behaviors (Duarte et al., 2021), in the extent of sibling cooperation (Rebar et al., 2020;Schrader et al., 2017), and in their larval morphology (Jarrett, Evans, et al., 2018).
At the 30th generation of experimental evolution, we used pooled whole-genome resequencing of these populations to document genetic variation at the molecular level when care was present and when it was prevented experimentally (Figure 1A).First, we determined the effect of care on within-population genetic variation (SNP diversity).Second, we identified the genetic loci that had diverged to the greatest extent following the removal of care by looking for regions of high genetic differentiation (F ST ) between experimental populations.

Breeding design and experimental evolution
We sampled DNA from experimental populations of N. vespilloides that had been evolving under different regimes of parental care and that were founded from a single genetically diverse population generated by interbreeding beetles from multiple wild populations across Cambridgeshire.These populations have been described in detail previously (Schrader et al., 2017) and comprise a total of four populations: two blocks (Block 1 and Block 2; separated by 1 week) containing two populations evolving with (FC POP ) or without parental care (NC POP ).Each replicate of each population originated from the same founding population and, therefore, is expected to be genetically identical.However, each block was bred 1 week apart and there could have been minor genetic and/or environmental fluctuations that could contribute to genetic variation within and between blocks (Barghi et al., 2019;Schlötterer, 2023).For the first 14 generations, when directional selection was high, an average of 34 pairs of unrelated beetles were bred at each generation (Schrader et al., 2017).Thereafter, populations were maintained with an average of 37 and 49 pairs at each generation for FC and NC populations, respectively, and were equivalently successful across the generations (Supplementary Table S1).On the 29th generation (as in every generation previously), we paired sexually mature males and females within each population.Each pair was placed in a separate breeding box with moist soil and a thawed carcass (10-12 g).We then placed each breeding box in a Figure 1.(A) Populations evolved in the presence (Full Care; FC) or in the absence of (No Care; NC) for 30 generations (two replicates per condition).Larvae were pooled for each replicate population (see Methods) for whole-genome sequencing (WGS).Distribution (top) and median (bottom) of (B) Pi (π) and (C) Watterson's theta (θ) across 1,000-bp nonoverlapping windows for FC and NC populations (error bars represent 95% bootstrapped confidence intervals).Block 1 (dashed line) and Block 2 (solid line) are plotted separately.cupboard and allowed parents to prepare the carcass and for the female to lay the clutch of eggs.For the NC POP , after 53 hr, both parents were removed from the nest just as had occurred for the prior 29 generations.Approximately 80 h after hatching we randomly selected 2-3 larvae from each family (15-18 families per population) for DNA extraction.

Larval tissue dissection, DNA extraction, and whole-genome sequencing
For each family, DNA from first-instar larvae were pooled and extracted using a modified version of the Qiagen DNEasy Mini Kit.Total DNA quality was checked using gel electrophoresis, and yield was quantified using a Qubit DNA Assay Kit (Thermo Fisher).Families were pooled in equimolar concentrations such that each individual was represented equally to generate four libraries: FC1, FC2, NC1, and NC2 with pool sizes of 41, 52, 52, and 59, respectively.Whole-genome resequencing libraries were constructed and sequenced (150-bp paired-end) at a depth of 100× using an Illumina Novaseq 6000 platform by Novogene (Hong Kong).

Intrapopulation genetic variation
We used π and Watterson's θ to measure levels of standing genetic variation within populations.Watterson's θ represents the expected number of segregating sites observed between a pair of homologous sequences sampled from a given population, whereas π is the average number of pairwise differences between all possible pairs of individuals in the sample.These measures were calculated for nonoverlapping 1,000-bp windows (for sites with coverage between 40 and 700 reads) across the genome using tools from Popoolation.We also computed genewise synonymous and nonsynonymous pi for CDS coordinates of all genes extracted from the reference annotation using Popoolation.
To allow comparisons to F ST windows we computed Tajima's D for 500-bp sliding windows with a 250-bp overlap.For all genetic diversity measures, we used non-parametric Kruskal-Wallis tests to test for differences between Full Care and No Care populations separately for each replicate block except in the case of Tajima's D. For normally distributed Tajima's D values, we used t-tests to test for differences between care conditions within each block.
Windows were filtered, so that statistics were based on windows that were covered across all replicates of both populations.

Genetic divergence between populations
To estimate population structure and demographic history, we extracted SNPs from the population sync file using the R package poolfstat (Gautier et al., 2022) using the core model of BayPass version 2.3 (Olazcuaga et al., 2020).Baypass uses allele frequencies to estimate a scaled covariance (Ω) matrix, which can be interpreted as the pairwise estimates of differentiation between the population.The Ω matrix was converted to a correlation matrix in R and visualized as a tree using the base R stats package.
To further measure the extent of genetic divergence between populations, we used Popoolation2 (Kofler, Pandey, et al., 2011) to calculate the pairwise fixation index (F ST ) for all combinations of population pairs across 500-bp sliding windows (250-bp overlap) across the genome.SNPs were called using sites with read counts between 40 and 700.Hierarchical clustering indicated that NC1 and NC2 were more closely related to their FC counterparts than to each other (Supplementary Figure S1).Moreover, inspection of F ST values across the genome indicated that the overall magnitude of differences between FC and NC differed between the blocks (Supplementary Figure S2).Therefore, to identify windows where evolving populations may have diverged consistently, over and above any variation within and between blocks, we computed F ST for each replicate line separately (i.e., FC1;NC1 and FC2;NC2).We then performed Fisher's exact tests for each of these windows to screen for significant allele frequency differences.We took the product of the −log(p) values for each block (FC1;NC1 x FC2;NC2) and selected the top 0.5% of values as regions of interest.In this way, we selected for loci which diverged consistently across the blocks, assuming that inconsistent divergence may reflect drift.Location of windows of interest were annotated using the reference genome and the intersect command in bedtools (Quinlan & Hall, 2010).A hit was considered only if the window intersected with the coordinates (either a gene or 5ʹ UTR) of the annotation by at least 1 bp.To characterize the extent of regulatory change, we took hits at annotated genes and further classified these genes into four possible categories (see Supplementary Table S3): (a) genes that encoded transcription factors (see Functional Annotation Methods for further information); (b) genes that encoded a gene involved in gene regulatory activities (gene expression-related; e.g., chromatin modifier, RNA polymerase, transcriptional cofactors, etc.); (c) long noncoding RNA; and (d) other protein-coding genes (those that did not fall into the first three regulatory function categories).
Using the same logic we used the auxiliary model in BayPass to identify candidate SNPs that were consistently associated with the loss of care across both blocks.Using the covariance structure among the population allele frequencies (Ω), the model explicitly accounts for the shared history of the populations, rendering the identification of SNPs potentially subjected to selection less sensitive to the confounding effect of demography (Gautier, 2015;Günther & Coop, 2013).Specifically, the model involves the introduction of a binary auxiliary variable to classify each locus as being associated or not with the loss of care.This allows the estimation of posterior inclusion probabilities (and Bayes factors [BF]) for each SNP while also accounting for multiple testing issues.For each SNP, the Bayes factor was expressed in deciban units (dB) via the transformation 10log10(BF).Significance was assessed based on the BF between models and SNP markers with strong evidence (BF > 20) were retained as potential candidates of interest (according to Jeffrey's rule) (Jeffreys, 1998).We then examined where these SNPs were located by looking for genes within 500 bp of the outlier SNP (using bedtools "window") making this comparable to our windowed approach.

Functional annotation
Functional enrichment analyses were conducted using the topGO R package version 2.38.1 (Alexa & Rahnenfuhrer, 2009) to identify overrepresentation of particular functional groups within the diverged genes in response to the removal of care, based on GO classifications using Fisher's exact test.GO terms were annotated to the N. vespilloides genome using the BLAST2GO (version 5.1.1)workflow to assign homologs to the Drosophila nonredundant protein databases (Gotz et al., 2008).To improve the GO term assignment, N. vespilloides genes were further annotated using a custom script that assigned GO terms from multiple well-annotated insect species (e.g., A. mellifera, B. terrestris, A. cephalotes, N. vitripennis, T. castaneum, and O. taurus) based on ortholog assignments obtained using Orthofinder (Emms & Kelly, 2019) using a custom pipeline (https://github.com/chriswyatt1/Goatee).To identify transcription factors we searched for the presence of known Pfam (Mistry et al., 2021) transcription factor domains in the protein sequences of the gene candidates of interest using Interproscan (Blum et al., 2021).Putative promoter regions (5ʹ UTRs) were classified as the 500-bp region upstream of each gene start coordinate (Quinlan & Hall, 2010).

Standing genetic variation between populations
First, we determined the effect of care on within-population standing genetic variation by measuring genetic diversity.We computed both Watterson's theta (θ) and Pi (π) statistics for each population across 1,000-bp nonoverlapping windows.Populations that evolved under Full Care (FC1 and FC2) had higher theta values than populations evolved under No Care (NC1 and NC2) (Table 1; Figure 1C; all p's ≤ .001).Similarly, there were higher Pi values in FC compared to NC, though this effect was not present in Block 1 (Table 1; Figure 1B).Together, these results suggest that FC populations maintained more SNP diversity compared with populations evolving under NC with some detectable variation between blocks (Table 1).We measured Tajima's D (500-bp overlapping windows) to further characterize the evolutionary forces shaping genetic diversity between populations.We show that genome-wide levels of Tajima's D are negative, with both replicates showing a significant reduction in Tajima's D in NC compared with FC (Table 2; Figure 3).

Genetic differences between populations
Next, we identified the genetic loci that had diverged to the greatest extent following the removal of care by looking for regions of high genetic differentiation (F ST ) between experimental populations (see Supplementary Figure S1 for population structure).We looked for changes over and above drift by looking for highly consistent divergence across the replicates using a 500-bp sliding window approach (see Methods; Figure 2A).Highly significant windows overlapped with both protein-coding and regulatory features of the genome, with 16% of windows being classified as a regulatory change in contrast to 47.9% of windows in protein-coding genes (Figure 2B).Using this approach, we identified 648 differentiated genes (Supplementary Table S3), with 144 of these windows uniquely intersecting the putative 5ʹ UTR regions of these genes only (Figure 2B; Supplementary Table S4).These genes were generally enriched for GO processes associated with morphogenesis, neural development, immunity, and hormone Table 1.Genetic diversity measures for each population evolving under Full Care (FC) and No Care (NC) for each block.Delta is the difference between FC and NC populations computed separately for each block (*p < .05,**p < .001).

Block
Evolving    S3) (B) cytochrome P450 6k1-like and (C) 5-hydroxytryptamine receptor-like.See Supplementary Figure S5 for replicate blocks plotted separately.All error bars represent 95% bootstrapped confidence intervals.signaling (Figure 2C; Supplementary Table S5).To test that our approach converged with other methods, we also identified SNP outliers using a Bayesian approach (see Methods).This method identified 3,086 outlier SNPs with consistent allele frequency differences between the NC and FC populations across both replicate blocks, which fell within 500 bp of 1,176 genes (Supplementary Table S6).These SNP outliers broadly converged on our windowed approach (220 genes; Supplementary Figure S3) with several key genes identified in both methods (Supplementary Table S7).
To test the hypothesis that genes selected in the No Care lost genetic variation, we examined the genewise ratio of nonsynonymous to synonymous π (π N /π S ; Table 1 and Supplementary Figure S4) and Tajima's D within 5 kb of divergent loci (Table 2; Figure 3).Both measures were reduced in No Care populations relative to the Full Care populations, in both replicate blocks, and this was true genome wide as well as for the genes identified in our divergence screens.

Discussion
We found that populations with parental care (Full Care) had greater levels of genetic variation, in the form of higher theta and pi diversity, than the populations where care was prevented (No Care).Previous work has suggested that social behavior contributes to genetic diversity mainly because of its effect on demography and particularly because of its influence on effective population size (Ne).Population genetic theory predicts that genetic diversity will increase with Ne and mutation rate (Charlesworth, 2009).Previous empirical work linking behavioral and life-history traits, such as reproductive strategy, fecundity, and body size with genetic diversity has suggested that these associations are ultimately mediated by changes in Ne (Bharti et al., 2023;Chak et al., 2022;Romiguier et al., 2014;Settepani et al., 2017).However, our results cannot be explained by demography because populations were maintained at similar population sizes with no differences in fecundity (Schrader, 2017) and no possibility of overlapping generations and/or changes in mating structure.Any small deviations in population size between care treatments were biased toward reducing genetic diversity in the Full Care populations-yet we found the opposite result.We suggest instead that the accumulation of genetic variation here is due directly to the effect of parental care in relaxing selection.The founding wild populations were inclined to provide care (Jarrett, Evans, et al., 2018) and likely had already accumulated high levels of standing genetic variation, which was swiftly lost when we exposed populations to selection in a No Care environment.
The majority of this accumulated genetic variation is likely to be either neutral or mildly deleterious, since the majority of new mutations generally fall into either of these two categories (Baer et al., 2007;Lynch, 2016).Indeed, we have previously demonstrated that inbreeding of these populations resulted in faster extinction of Full Care compared with No Care populations, further suggesting that, at least some, of the variation accumulated in the presence of care was deleterious (Pascoal et al., 2023).Although we measured only SNP variation here, genetic variants (e.g., insertions/deletions, transpositions) that arise through different types of mutation or recombination could also, in theory, be maintained in the population by parental care.Whether care favors particular types of mutants remains to be tested in future studies.In contrast, the harsher No Care environment imposed strong directional selection resulting in rapid adaptation (Schrader et al., 2015(Schrader et al., , 2017) ) and reducing levels of standing genetic variation.We identified genetic divergence at number of loci, which were also associated with the loss of nonsynonymous pi (lower π N /π S ) and reductions in Tajima's D, a pattern that was similar to the genome-wide differences in genetic diversity.Again, this is consistent with the interpretation that No Care populations experienced strong directional selection, while Full Care populations harbored more potentially deleterious mutations.
Here, we follow convention in assuming that loci that diverge consistently are likely to represent adaptive genomic change, whereas inconsistent responses to the parental care treatment (No Care vs.Full Care) are due to drift.Yet alternative explanations for inconsistent responses are also possible, and this might be particularly true for traits under social selection (as opposed to abiotic selection pressures).Inconsistent patterns of genetic change across experimental blocks attributed to drift might instead reflect idiosyncratic or opportunistic responses to selection that arise through subtle variation in founding populations (Barghi et al., 2019;Brennan et al., 2022;Schlötterer, 2023).This is not surprising given that polygenic traits can be genetically redundant and adaptation can arise through multiple intersecting pathways and unique combinations of alleles within a population (Barghi et al., 2019;Láruson et al., 2020).Moreover, the magnitude of these inconsistent differences, due to either drift or selection, might have been intensified by the selection regime imposed by the social environment, depending on whether it relaxed selection or imposed directional selection, for example, or whether the strength of selection was modulated by genes of social partners, parents, or siblings (Drown & Wade, 2014;Linksvayer & Wade, 2009).Such effects could explain variation within and between replicate populations that accumulates over time.Although we cannot distinguish idiosyncratic adaptive change from drift with our data currently, future work using a high number of replicated populations measured across several generations could provide key insights into these evolutionary dynamics (Barghi et al., 2019).
Our data suggest that No Care populations diverge from Full Care populations at loci that could promote immunity, metabolic, and behavioral stress resilience in the absence of care.The loss of care in N. vespilloides is likely to be associated with greater levels of environmental stress during development and heightened exposure to pathogens from the carrion resource (Mashoodh et al., 2021;Rozen et al., 2008).We have previously shown that adaptation to a No Care environment is associated with gene expression signatures that show blunted stress responses and compensatory expression in metabolic and developmental pathways (Mashoodh et al., 2021).Not surprisingly, many of the genetic differences between the populations are in upstream regions and/or genes that encode for transcription factors, chromatin modifiers, and other genes that modify transcription, suggesting that change in regulatory function through varied mechanisms is a key component of adaptation to the loss of care.This is likely to be an underestimate of the extent of regulatory change, as we have yet to characterize the regulatory landscape of N. vespilloides and windows without an annotated overlap could be in distal promoter and/or enhancer regions.Nevertheless, differences in regulatory functions could shape levels of gene expression of other genes, further buffering against stress in the absence of parental care (Mashoodh et al., 2021).In this way, the signatures of adaptation to the loss of parental care are not much different to adaptive genetic responses to other abiotic stressors in the broader environment.Indeed, a key feature of stress adaptation across species is that it involves changes in gene regulatory pathways and this is true from bacteria to plants and animals (Barghi et al., 2019;Bhargava & Sawant, 2013;Conrad et al., 2010;De Nadal et al., 2011).Although we cannot identify a single gene or master regulator within the regulatory changes, these data do identify candidate regulatory genes that might play key roles in conferring resilience to the loss of care, and to environmental stressors more broadly.
Delving more deeply into loci at which we detected the greatest differences, we found that adaptation to the loss of care involved changes in several genes associated with immune function (e.g., CD109 antigen and lysozyme c-1), which could help cope with the increased exposure to the bacterial pathogens of the carcass nest experienced by No Care larvae.Previous work on N. vespilloides has shown that lysozyme expression is particularly heightened in parents immediately after the larvae hatch (Cotter & Kilner, 2010;Palmer et al., 2016) and that it is likely to be particularly important for eliminating pathogenic Gammaproteobacteria (Duarte et al., 2018).Relatedly, we found divergence at a number of cytochrome P450 genes (4ac1, 4ac2, 4c1, 4g15, 9e3, and 6k1-like; Supplementary Table S3), which are known to be involved in the metabolism of endogenous compounds as well as exogenous toxins and disease vectors, and which might also participate in defensive responses (Nauen et al., 2022).
Cytochrome genes also appear to play a role in mediating the response to social density in Drosophila.Expression changes at the Cyp4, Cyp6, and Cyp9 gene families can be induced by manipulating social density in Drosophila and deletions of the Cyp6a20 gene have been associated with higher levels of aggression and reduced sociality (L.Wang et al., 2008).This is particularly interesting given that we have previously shown that larvae from the No Care populations evolved to show greater levels of sibling cooperation than larvae from the Full Care populations (Rebar et al., 2020;Schrader et al., 2018).Cytochrome P450 families tend to share similar functional domains, and therefore, it is possible that changes in these genes could have effects on social behavior via their actions on multiple hormonal systems (e.g., pheromones, ecdysone) (Iga & Kataoka, 2012;Nauen et al., 2022).Furthermore, P450 genes are also intertwined with juvenile hormone pathways, which are known to be involved in multiple facets of behavioral and morphological development (Flatt et al., 2005).This could include adaptations that aid in locating and facilitating the use of the carrion breeding resource, such as the increase in relative mandible size and reduced arrival time at the carcass that we also detected in No Care larvae (Jarrett, Evans, et al., 2018;Jarrett, Rebar, et al., 2018).This interpretation is additionally consistent with most of the genetic hits belonging to cell signaling and biosynthetic pathways, which fall into GO categories associated with morphological, brain, and olfactory development.
We also found changes in neuropeptides that are also involved in metabolic, homeostatic, and feeding pathways (e.g., orexin, 5-hydroxytryptamine, and cholecystokinin receptors), raising the possibility that these could represent adaptations in larvae for feeding and extracting nutrients from the carcass resource in the absence of parents (Nässel & Zandawala, 2019).Previous work has shown that genes ancestrally associated with metabolic, homeostatic, and feeding pathways can be co-opted to serve new social functions (Potticary et al., 2022(Potticary et al., , 2023)).For example, the oxytocin/vasopressin system is commonly associated with the expression of parental care, pair bonding, and other affiliative social behaviors in mammals (Froemke & Young, 2021), but has an ancestral role associated with promoting water balance (Koto et al., 2019).A recent study in the closely related burying beetle, N. orbicollis, suggested that the expression of inotocin (the insect homolog of oxytocin/vasopressin) was correlated with the transition to parenting, an effect that was more pronounced in males than in females (Potticary et al., 2022).Takeout is another gene that, despite being typically associated with feeding and circadian rhythms, has been shown to be highly expressed while parenting in burying beetles and may be involved in the transition from infanticide to larval care (Moss et al., 2022;Potticary et al., 2023;Sarov-Blat et al., 2000).We found divergence in both an oxytocin receptor and a takeout homolog (Supplementary Table S3), which could explain how male parental care eventually decayed in the No Care lines (Bladon et al., 2023).Finally, angiotensin converting enzyme was one of the most diverged genes in our analyses.This gene has varied roles from conferring immunity to regulating neuropeptide signaling (Isaac et al., 1999).However, it also appears to be strongly expressed in insect reproductive tissues and may, therefore, play a role in adaptations in mating and fecundity between the populations (Schrader, 2017).
While we show here that parental care contributes to genetic variation through its effect on selection, it is possible that the incidence of mutation is itself reduced by the loss of parental care.Mutation rates have a strong genetic basis and can vary between individuals and among populations (Baer et al., 2005;Y. Wang et al., 2023).Given that the loss of care is a major developmental stressor, and that stress has been shown to induce mutations, adaptation to the loss of care could involve genetic mechanisms that dampen and/or buffer the consequences of new mutations that arise (Baer et al., 2007;Snell-Rood et al., 2016).Consistent with this hypothesis, we found high levels of genetic differentiation among genes involved in DNA replication and repair (e.g., Artemis and the PAXIP1 interacting protein; Supplementary Table S3) (Kurosawa & Adachi, 2010;Muñoz & Rouse, 2009).These genes, as part of their role in stress regulation, could facilitate efficient DNA repair, purging new mutations and shaping the subsequent mutation load of a population.The observation that multiple transfer RNAs (tRNAs) show divergence is particularly interesting given that variation in tRNAs has been associated with increased mutation loads via transcription-assisted mutagenesis (Thornlow et al., 2018).In other words, both genetic and phenotypic adaptations within each population could favor an optimal mutation-selection balance, resulting in different levels of standing genetic variation based on levels of care experienced within each population.Further functional characterization of these changes would help clarify if parental care facilitates the evolution of the mutation rate, potentially providing another mechanism for divergence in mutation-selection balance among populations (Baer et al., 2007;Lynch et al., 2016).
In short, we have shown that parental care allows genetic variation to accumulate by relaxing selection.When care is lost, a number of genetic changes quickly follow, which may be adaptive and which result in the loss of standing genetic variation.Better functional characterization of these gene targets and regulatory regions is now required to understand the genetic causes and functional consequences of the differences we have found between populations that are, and are not, exposed to posthatching care.These are key areas of future work that will help explain whether the maintenance of standing genetic variation under parental care is likely to help or hinder adaptation in a rapidly changing world.

Figure 2 .
Figure 2. (A) Product of -log(p) values for Fishers exact tests between Full Care (FC) and No Care (NC) populations of each block (i.e., FC1;NC1 x FC2;NC2) for each 500-bp sliding window (250-bp overlap) sorted by position.Dashed red line indicates 99.5th percentile.(B) Percent of windows overlapping with genomic features (ncRNA = noncoding RNA; 5ʹ UTR is defined as 500-bp upstream of gene start position) and the number of genes that correspond to each gene category (see Methods).(C) Representative enriched GO terms (biological processes) for the most diverged genes between FC and NC populations.

Table 2 .
Tajima's D (mean of 500-bp sliding windows) for populations evolving under Full Care (FC) and No Care (NC) for each block.Statistics are presented for all windows across the genome (genome-wide) as well as for windows that overlapped with diverged genes (± 5 kb).Delta is difference between Tajima's D between means of FC and NC populations (**p < .001).