Heredity

Demographic inference based on sequentially Markovian coalescent (SMC) model is widely used for disentangling the evolutionary history of species and populations. However, the genetic effects of self-fertilization, including the expected reduction of population-level genetic diversity and hence effective population size (Ne), are often ignored in demographic inference, despite the fact that selfing is a common reproductive mode in both invertebrates and plants. To investigate how self-fertilization could bias demographic inference, we apply eSMC2, a selfing-aware demographic inference method, to Hawaiian populations of Caenorhabditis elegans, a predominantly self-fertilizing nematode. Using an expanded dataset that contains more Hawaiian strains than in previous studies, we recovered previously reported genetic groups, including one largely consisting of strains from the Big Island. Population substructure is also detected among these Big Island subpopulations, highlighting the importance of identifying subgroups to avoid misinterpretations of demography. For different Big Island subpopulations, we estimate high selfing rates compatible with previous studies and infer consistent population declines. Population sizes inferred using eSMC2 are smaller than those inferred using PSMC', which does not account for selfing. However, rescaling PSMC' parameters to account for self-fertilization produces a similar demographic trajectory as eSMC2 inference. Our study thus demonstrates how one can account for effects of selfing when using SMC-based methods. It also indicates that the Hawaiian group shows population declines similar to those previously found outside Hawaii.

Arctic ecosystems are undergoing rapid climate warming, providing an opportunity to examine molecular evolution under environmental stress. Here, we used whole-genome resequencing to assess population structure, genetic diversity, and recent demographic history in the snowy owl (Bubo scandiacus), a highly mobile Arctic apex predator classified as Vulnerable by the IUCN. We analysed genome-wide data from individuals sampled across the species' circumpolar range using a chromosome-level reference genome. Population genomic analyses revealed an absence of detectable population structure, indicating extensive gene flow and near-panmixia across continents. Genome-wide nucleotide diversity was low ({pi}= 2.47 x 10^(-4)), while runs of homozygosity were dominated by short tracts, indicating limited standing genetic variation despite little evidence for recent close inbreeding. Our demographic inference revealed a sustained decline in effective population size over the past ~1,600 years, a pronounced reduction coinciding with the Medieval Warm Period followed by a partial recovery during the subsequent cooler interval. This direct, measurable effect of past climate warming on effective population size raises concerns for the long-term survival of the snowy owl. Our results demonstrate that high connectivity and low levels of recent inbreeding have not necessarily buffered this species against climate-driven demographic erosion of genetic diversity. This study highlights the value of using whole-genome population genomics to inform genomic monitoring and conservation decision-making in rapidly changing environments.

Amongst measures of linkage disequilibrium, the correlation measure r2, due originally to Sewall Wright, has several advantages. Although basically a haploid statistic, it can readily be estimated using covariance and correlation from diploid data. Its population expectation can also be found from probability methods. We show that the probability of linked identity by descent, LIBD, directly estimates the expectation for r2 for Wright-Fisher models. The simplest such model for two loci is the haploid model, where each haplotype is chosen independently from haplotypes of the previous generation. This gives an equilibrium between recombination (recombination frequency c) and drift (population size N) where the expected value of r2 is approximately 1/(1+4Nc). For diploid models, the expectation is similar for low values of c, but contains a term in c2 and rises to double the haploid value for the limiting value of c = 0.5. An increase in r2 is also expected for the model of monogamous mating, as opposed to random mating for each offspring. The expectation for r2 is unaffected by the use of diploid rather than haploid data, an exception being the case of monogamous mating. Computer simulation shows close agreement of r2 with probability expectations, with minor discrepancies observed for unlinked loci. Using the probabilistic definition, r2 can be intuitively understood as a parameter shaping the effective sample size in genome-wide association studies.

Demographic history can profoundly shape genomic variation by altering levels of inbreeding, the strength of genetic drift, and the efficacy of purifying selection, thereby influencing both the accumulation and expression of deleterious variation. Yet the extent to which these genomic consequences translate into organismal performance remains poorly understood. Here, we investigated how contrasting demographic histories shape genomic diversity, inbreeding, genetic load, and performance across mainland and insular populations of the Italian wall lizard, Podarcis siculus. Using whole-genome sequence data from 263 wild-caught individuals and field-based performance measurements from 134 lizards, we characterised population structure and demographic history, quantified genomic diversity, runs of homozygosity, and deleterious variation, and tested associations between genomic inbreeding and performance traits. Island populations showed reduced genetic diversity, increased autozygosity, and elevated genetic load relative to mainland populations, consistent with stronger genetic drift and reduced efficacy of purifying selection. Patterns of deleterious variation were consistent with purging of strongly, but not mildly, deleterious variants. Greater genomic autozygosity was associated with reduced bite force in island lizards, a pattern consistent with inbreeding depression, whereas sprint speed, sperm quality, immunocompetence, and spatial cognition showed no significant association with inbreeding. Genome-wide mapping of runs of homozygosity further identified multiple genomic regions associated with reduced bite performance. Together, these results show how demographic history shapes genomic diversity, homozygosity, and deleterious variation in wild vertebrate populations, while revealing that phenotypic consequences can be trait-specific, a pattern that may reflect differences in genetic architecture and the extent of environmental influence on trait variation.

Extreme climatic events (ECEs) could potentially induce episodes of rapid adaptive evolution, altering the evolutionary trajectories and viability of wild populations experiencing increasingly stochastic environments. Adaptive potential stems substantially from additive genetic variance (V) in survival. However, the degree to which ECEs impact this variance, and the resulting evolutionary implications, have not been fully conceptualised or quantified in nature. These aims are highly challenging because underlying genetic and environmental variances have fundamentally non-linear effects on the dichotomous phenotypic outcome of survival versus death, while substantial data are required for estimation. We present a conceptual and analytical framework that formalises biological hypotheses for how V depends on environmental severity, how such dependences are reshaped by non-linear genotype-environment-phenotype relationships and demographic structure, and the implications for evolutionary responses. We then use 25 years of survival and pedigree data from European shags (Gulosus aristotelis) to decompose variance in adult survival across diverse environmental conditions including ECEs. Using a capture-recapture animal model, we demonstrate non-negligible V underlying adult survival across all years combined, implying notable adaptive potential. Yet, partitioning by environmental severity revealed substantial V in good years with high survival but negligible V in bad years with ECE-induced low survival, implying a genotype-by-environment interaction that constrains adaptive potential under severe conditions. However, scale transformations showed how the non-linear genotype-environment-phenotype relationship can distort or even reverse such patterns, by reducing the extent to which V is phenotypically expressed in good conditions with high mean survival and amplifying it in more severe conditions with lower mean survival. Evolutionary responses to severe conditions could consequently be non-negligible even given small underlying V. We thereby demonstrate how ECEs can reshape adaptive potential, highlighting the need to integrate impacts of environmental severity on genetic variation in fitness components into predictions of population resilience under climate change.

Understanding the demographic processes shaping population connectivity is essential for conservation. Yet disentangling the effects of dispersal and local effective density remains challenging because traditional genetic approaches integrate signals over evolutionary timescales. Identity-by-descent (IBD) segments provide opportunities to reconstruct recent demographic connectivity over timescales relevant to contemporary conservation challenges, but their reliability under realistic ecological conditions remains poorly evaluated. Here, we assess the ability of two promising spatiall IBD-based methods, IBD-Analysis and MAPS, to estimate effective density and dispersal variance. Using a simulation-inference framework combining spatially explicit simulations, IBD segment extraction, and demographic inference, we evaluated both methods across a wide range of demographic, ecological and methodological scenarios. Both approaches reliably recovered demographic parameters under homogeneous landscapes, with IBD-Analysis providing accurate estimates of global effective density and dispersal across a broad range of conditions. Spatial heterogeneity emerged as the main driver of inference performance. When demographic parameters varied across space, IBD-Analysis integrated local variation into biased global estimates, whereas MAPS captured spatial patterns but showed substantial variability and generated spurious spatial structure in homogeneous landscapes. Methodological factors had comparatively weaker effects than ecological complexity, highlighting the importance of accounting for landscape structure when applying IBD-based demographic inference. These results highlight the potential of IBD-based approaches to move beyond descriptive patterns of genetic connectivity toward mechanistic inference of demographic processes. By distinguishing the contributions of dispersal and density, these methods can improve conservation decisions by linking genomic patterns to key demographic parameters relevant for conservation.

Polyandry, multiple mating by females, both benefits females by increases genetic variation within families and provides an opportunity for sexual selection through sperm competition and cryptic female choice. Estimating polyandry in species with large clutches is challenging because genotyping all progeny individually can be impractical. We present a likelihood-based framework for estimating the number of sires and their relative contributions from the maternal genotype and pooled sequencing of entire progeny pools. The method evaluates the probability of observed pooled offspring read counts under alternative paternal genotype mixtures, accounting for Mendelian segregation and sampling variation from finite sequencing depth. We implement this framework for both biallelic SNPs and short-read mini-haplotype markers constructed from closely linked SNPs, which provide multi-allelic information comparable to traditional microsatellite markers. Simulations show that this approach accurately recovers sire number and sire contributions across a range of conditions. We then applied the method to 275 wild-caught Drosophila melanogaster females and their offspring from a Kansas population using 3,053 biallelic SNPs and 239 four-allele mini-haplotype markers. Estimates from the two marker sets were highly concordant, and combined analyses indicate that 57% of families were sired by a single male. These results demonstrate that pooled progeny sequencing provides an efficient genomic approach for estimating polyandry in species with large families.

Habitat fragmentation and destruction in old-growth forests represents a primary threat to specialist species, frequently resulting in the erosion of genetic diversity. The Hermit Beetle, Osmoderma barnabita, is a protected relic of primeval forests and an obligate saproxylic specialist whose larval development depends exclusively on the wood-mould of hollow trees. This study provides the first high-resolution, genome-wide assessment of O. barnabita across seven Polish forest complexes using double-digestion Restriction-site Associated DNA sequencing data from 49 individuals. Following rigorous quality control and reference bias assessment, we evaluate contemporary population structure. To distinguish between evolutionary processes, we isolated 2,920 putatively adaptive outliers from the neutral background, evaluating both neutral and adaptive genomic signals. The analyses revealed four distinct genetic clusters, confirming that the Polish metapopulation is severely subdivided. Crucially, the westernmost population in the Lower Silesian Forest emerged as the most divergent evolutionary unit (FST up to 0.360), suggesting potential hybrid origin or long-term isolation. We observed a striking linear correlation (R2=0.986) between neutral diversity and adaptive potential, indicating that genetic drift in small and isolated forest fragments is actively depleting the variation required for long-term persistence of this species. Our findings demonstrate significant gene flow restriction across all sites, suggesting that current conservation frameworks fail to encompass these highly differentiated lineages. To mitigate this genetic vulnerability, we recommend reclassifying these clusters as independent Management Units and Evolutionary Significant Units. Maintaining forest continuity and protecting local micro-habitats are essential to prevent further isolation and preserve the genomic integrity of this saproxylic flagship species.

Demographic stochasticity is generally expected to contribute little to population variability when census abundance is large, because random demographic events are increasingly averaged across many individuals. This expectation, however, depends on how evenly individuals contribute to reproduction. In many natural populations, reproductive success is highly skewed, so that a large census population can coexist with a much smaller effective number of breeders (Nb). Although Nb is commonly used as a genetic measure of contemporary effective size, it also summarizes variation in reproductive contribution and may therefore contain information about demographic variability. Here, we develop a female-based age-structured theoretical framework linking Nb to temporal variability in mature population abundance, with N denoting the census number of mature females and Nb the effective number of female breeders. Recruitment variability arises from non-Poisson reproductive variation and environmental stochasticity, and its transmission to mature abundance is described by a life-history coefficient, {gamma}, determined by survival, maturation, and the contribution of overlapping cohorts. When the Poisson component of recruitment variance is negligible and Nb/N is sufficiently small, the coefficient of variation in mature abundance, CV[N], can be approximated from Nb, environmental variability, and {gamma}, with little additional information from census abundance. This approximation does not require specification of the functional form of the relationship between adult abundance and recruitment. Stochastic age-structured simulations closely reproduced the analytical prediction, even when additional demographic stochasticity in survival was included. Substantial fluctuations in mature abundance occurred at census sizes on the order of 105 when Nb was small, whereas overlapping cohorts buffered recruitment fluctuations. Because Nb can be estimated genetically from samples of a single cohort, the framework provides a direct link between genetic monitoring and population dynamics. More broadly, these results show that large census abundance alone does not guarantee weak demographic variability and highlight the potential of genetic data to inform expected population fluctuations when detailed demographic information is limited.

Social genetic effects can substantially influence phenotypes in animals and plants and the study of them is of interest. Group housing in animals for example, in pigs provide an ideal framework for studying interactions among individuals. In animals, the genetic of social effects has been studied by simplifying the underlying interactions for example, by assuming that all the individuals in the group contribute equally to the social environment. This homogeneous behaviour is not observed in real conditions. Here, we address this issue showing that identification of interacting individuals is feasible in pig behaviour during growth, and we propose a model that leads to more accurate modelling of social genetics. Our approach infers the strength of interaction via time interval between visits of animals to automatic feeders. This measure of strength was then used in a reaction norm model to evaluate the genetic determinism of social effects. We also build a framework to transform variance components estimated for a trait continuously measured, such as feed intake per visit, into variance of total or aggregated trait e.g. during the growth period. Including social effects improves the fitting of the models. Visit-based social interaction was highest at short time intervals, decreasing rapidly and becoming negligible after 80 seconds. Variance components for the aggregated trait the total feed intake support the contribution of the social interactions. Our approach allows more precise understanding of social effects through individual feeding events and a measure of interaction strength, and further translates these into in an aggregated phenotype of interest.

Self-fertilization results in reduction of heterozygosity over time, accelerated by small population sizes, restricted ranges, and fragmented habitats. Yet some species maintain heterozygosity despite these factors. One is Corallorhiza bentleyi, a threatened, mycoheterotrophic, self-pollinating orchid restricted to Virginia and West Virginia, USA. Reduced representation sequencing was employed across five of the ~12 known localities of C. bentleyi, including representatives of sister species C. involuta from Mexico, to quantify patterns of population structure and genomic diversity. Surveys of C. bentleyi revealed a predominance of the cleistogamous (closed-flowered) growth form. Differentiation was observed between the two species but with weak population structure in C. bentleyi. Surprisingly, excess heterozygosity and negative inbreeding were recovered in both species, counter to expectations. Analyses provided evidence against bioinformatic artifacts (collapsed paralogs), small-sample bias (Robertson effect), and accumulation of mutations along chromosomes due to low recombination associated with selfing (Meselson effect) as explanations for excess heterozygosity. Demographic analyses favored a model of post-glacial divergence among C. bentleyi and C. involuta lacking bottlenecks. Finally, simulations under alternative selection scenarios favored models of overdominance vs. those of negative frequency dependent selection. It was concluded that overdominance, namely selection against selfed homozygotes, is the most likely explanation for elevated heterozygosity in C. bentleyi. Conservation concerns are discussed in light of these findings and the threats of extinction for one of the most vulnerable orchid species in North America.

Chromosomal inversions are proposed to aid local adaptation to divergent habitats when there is gene flow by creating regions of low recombination that can maintain tight linkage between multiple locally adaptive alleles. However, key theoretical predictions about the role of chromosomal inversions in local adaptation remain largely untested in empirical systems. Here, we test the prediction that inversions will concentrate relatively small genetic effects across multiple ecologically relevant traits by exploring inversion-phenotype associations of three inversions in threespine stickleback (Gasterosteus aculeatus). Using complementary approaches of association mapping in wild-caught individuals and genetic mapping in lab-raised individuals from two divergent genomic backgrounds, we found that an inversion on chromosome XXI influences multiple ecologically relevant traits. The inversion has small-to-moderate effects on these traits, as expected for inversions involved in polygenic adaptation. Further consistent with theoretical predictions, inversion heterozygotes show intermediate phenotypes. While most inversion-phenotype associations are consistent across different genomic backgrounds and experimental approaches, some inversion-phenotype associations were specific to a single dataset. These results indicate that the phenotypic effects of inversions are shaped by both genomic and environmental context, which could facilitate population-specific adaptation. Together, these results support the role of the chromosome XXI inversion in ecological divergence in threespine stickleback. More broadly, our results suggest that inversions can contribute to local adaptation by preserving linkage between small-effect adaptive alleles despite extensive gene flow.

Following a severe 2018 drought affecting much of Europe, European beech (Fagus sylvatica) exhibited patchy canopy dieback. To investigate the drivers of this mortality, we used a paired case-control design across Swiss forests, integrating phenotypic assessments and genomic analyses. Phenotypically, crown transparency and discoloration were used to classify tree damage state, while several other characteristics were measured to evaluate whether paired, neighboring damaged and undamaged trees were otherwise comparable, and to identify physiological stress signatures. Elevated leaf calcium content and red edge reflectance in dried leaves emerged as robust indicators of physiological stress. Genomic association models using reference-based SNP markers and reference-free k-mers were heavily confounded by micro-environmental variance and local population stratification even within our geographically limited study region. Our cross-validated genomic outliers showed no spatial overlap with predictive drought-associated loci previously identified in an independent German cohort. These findings indicate that drought survival in European beech is not governed by a universal set of major-effect alleles. Instead, tree performance appears highly context-dependent, likely shaped by micro-environmental heterogeneity and phenotypic condition. However, our results do not exclude polygenic drought responses involving multiple small-effect loci.

Populations adapting to heterogeneous environments often harbor many alleles that covary with environmental gradients, yet how so many adaptive variants can be maintained and show signatures of local adaptation remains unclear. Epistatic non-additive interactions among alleles may concentrate fitness effects at specific allelic combinations to shape complex traits and facilitate adaptation. We evaluate compositional epistasis among three intronic SNPs in the couch potato (cpo) gene, each independently clinal along the North American latitudinal gradient, using recombinant outbred populations fixed for distinct three-locus combinations. Flies carrying all three temperate alleles (cpoTTA) or all three tropical alleles (cpoCGT) show non-additive increases in traits associated with fitness in their respective environments: dormancy, starvation resistance, and lipid storage in temperate flies and fecundity in tropical flies. Intermediate genotype combinations generally show weaker and negative deviations from additive predictions. Temperate genotypes show reduced cpo expression and independently create predicted transcription factor binding sites, suggesting a regulatory basis for the observed phenotypic effects. Consistent with selection shaping haplotype frequencies against a backdrop of ongoing recombination, the temperate and tropical genotypes are enriched above random expectation across the latitudinal cline, with directionally consistent but non-significant seasonal patterns. Our findings demonstrate that non-additive interactions among non-coding variants within a single gene can produce large, pleiotropic fitness effects and shape genotype distributions in natural populations, providing a mechanism by which a small number of regulatory variants can underlie complex climatic adaptation.

Environmental exposures can induce epigenetic modifications that persist across generations, potentially contributing to the transmission of environmentally-induced phenotypes. However, the extent to which such persistence of molecular changes are maintained independently of genetic mechanisms remains unclear. In this study, we investigated the evolution and genetic determinism of DNA methylation (DNAm) variation across three successive generations of Japanese quails (G0 to G2) belonging to two epilines, defined according to whether their female ancestor had received genistein supplementation (epi+) or not (epi-). Using reduced representation bisulfite sequencing (RRBS), we characterised CpG methylation patterns and identified differentially methylated cytosines (DMCs) and regions (DMRs) between control and genistein-supplemented epilines. Out of 112,745 CpG sites analysed, no DMCs or DMRs were detected between epilines in the first generation following exposure (G0), whereas 621 DMCs (36 DMRs) were identified in G1, and 1,381 DMCs (101 DMRs) in G2. Despite this progressive increase of differential methylation sites, mean differences in methylation rate between epilines were globally small. Similarly, genome-wide genetic differentiation between epilines was very limited. Overall, SNP-based heritability of DNA methylation was 0.19 indicating that genetic variation accounted for only a small part of the global variation of DNAm levels. By contrast, DNAm at differentially methylated CpG sites was extensively under genetic control, with a mean SNP-based heritability > 0.5, and meQTL analysis further identified significant associations between SNPs and methylation levels. Together, these results indicate that local DNAm differentiation in quail is under substantial genetic regulation, while the limited genetic differentiation between groups suggests that genetic variation alone can not fully explain the progressive accumulation of methylation differences following an ancestral exposure to genistein. These results support DNAm as a plausible molecular candidate mechanism in the multigenerational response to environmental exposure and highlight the complex interplay between genetic and epigenetic regulations in the inheritance of phenotypes.

Populations that have detectable additive genetic variance for fitness (VA(W)) have the capacity for increase in their mean fitness through natural selection. The extent to which the magnitude of VA(W) in wild populations is predictive of the rate of adaptation is not well known, because few empirical field-based studies have been designed to address this question. Here, we report a two-year field experiment with pedigreed seeds embedded in four natural populations of the annual wildflower Baby Blue Eyes (Nemophila menziesii H. & A.), a species native to western North America. We first estimated VA(W) in a parental generation and then evaluated whether that estimate accurately predicted the adaptive response observed in the subsequent offspring generation grown under natural field conditions. We achieved this by directly comparing the mean fitnesses of the parental generation and their progeny grown under the same environmental conditions. We detected VA(W) and evidence of adaptation in all populations despite inter-annual differences in environmental conditions that directly affected fitness. Genotype-by-year interactions in fitness indicate that the fitness ranking of families differed between years, implying that the selective environment shifted over time. Additionally, although the environment in the second year reduced population mean fitness in three sites, the offspring generation maintained mean fitness at or above replacement, suggesting that adaptive evolution may have buffered the effects of environmental change. Overall, our results suggest that the four populations harbored sufficient VA(W) for adaptation through natural selection to the environmental conditions observed during the two years of this study.

Early embryos must mount stress defenses without compromising the developmental programs that drive morphogenesis, but whether genetic background alters this balance remains unclear. We compared cellularizing Drosophila melanogaster embryos from inbred and outcrossed backgrounds under control conditions and after acute or chronic heat-stress. Both backgrounds developed successfully and showed modest transcriptomic differences at 18{degrees}C. Heat stress at 32{degrees}C, however, sharply separated their phenotypic and transcriptional responses. Inbred embryos showed poor larval hatching after acute and chronic heat, and this deficit persisted through adult eclosion. Their transcriptomes exhibited a greater stress-response burden, with enrichment of proteostasis, DNA repair, mitochondrial maintenance, RNA processing, and metabolic remodeling programs. Outcrossed embryos retained substantially greater developmental survival and higher abundance of transcripts associated with development, adhesion, epithelial organization, signaling, transport, sensory responses, and morphogenesis. Background-dependent divergence was greatest during acute heat. Under chronic heat, background-dependent differences were fewer, with most overlapping those seen during acute heat and a smaller set detected only under chronic heat. Transcript-origin analysis further showed that acute outcrossed-biased genes were predominantly maternally derived and showed no enrichment for purely zygotic transcripts, whereas chronic outcrossed-biased genes contained a larger zygotic component. In summary, the less viable inbred background mounted the broader molecular stress signature, whereas the more viable outcrossed background retained greater abundance of development-associated transcripts. These findings identify genetic background as a major determinant of how embryonic heat responses are distributed between molecular stress defense and developmental gene expression.

Motivation: Optimum contribution selection (OCS) maximizes genetic gain under a constraint on relatedness. Genomic relationship matrices are dense, so the quadratic term costs O(n^3) operations, and existing OCS software assumes pedigree sparsity. Results: We replace the relationship matrix by its leading eigenpairs, computed by randomized SVD, and solve OCS in factorized form. On Norway spruce population (n=5,525), ranks 10--30 kept the exact objective within 1% of the optimum at 110--330-fold speedups; on the QTL-MAS 2010 (n=3,226), ranks up to 100 were needed (24-211 fold). The required rank rises as the diversity constraint loosens. Gain alone is a misleading accuracy check, since coancestry exceeded the optimum by up to 28%, so we choose the rank on the exact objective, which can be checked without a dense solve. The speedup comes from computing only the leading eigenpairs. We extend the factorization to robust OCS (18-fold). Availability and implementation: Julia implementation, and examples are available at https://github.com/jonhar97/fastocs under the MIT license (archived at Zenodo, DOI 10.5281/zenodo.23214804).

Signals of positive selection around the SLC30A9 gene have been reported in human populations outside Africa. Selection likely acted on a highly differentiated single-nucleotide polymorphism, rs1047626, leading to a non-synonymous substitution in the encoded zinc transporter. Because of the striking similarity between the putatively selected SLC30A9 haplotype observed in several current human populations and the Denisovan individual, previous work has proposed adaptive introgression. Yet alternative explanations, including ancient human variation, and the precise archaic source -Neanderthal or Denisovan- remained unresolved. Considering the potentially complex evolution of SLC30A9, we applied Approximate Bayesian Computation (ABC) algorithms coupled to machine learning to investigate the most plausible evolutionary origin of this substitution. After modelling different evolutionary scenarios with forward-in-time simulations, our results highlight that the most probable scenario is a Denisovan origin of the rs1047626 polymorphism. However, the allele likely introgressed into Neanderthals first and was then passed into non-African modern humans. Moreover, the derived allele frequency for rs1047626 across several African populations is consistent with back-to-Africa migrations. Finally, our ABC analyses indicate strong positive selection in East Asian populations and other out-of-Africa populations, whereas in Neanderthal populations, the selection coefficient was probably neutral or slightly deleterious.

The phenotype of a focal individual can be affected not only by the individuals own genes (direct genetic effects, DGE) but also by genes of other individuals it interacts with (indirect genetic effects, IGE). Better understanding IGE can help better understand the impact and mechanisms of social effects. Here, we investigated, in laboratory mice, the modulation of IGE by genetic background of the focal individual (social epistasis) and by sex. We focused on IGE arising from Epha4, a gene that has been associated with the individuals own stress-related and social behaviours (i.e. DGE) and the stress-coping behaviour and wound healing of the other mice in the cage (i.e. IGE). Among four different genetic backgrounds for focal mice, significant IGE from Epha4 were detected in CBA females only. A formal model selection analysis was consistent with the modulation of IGE by the similarity of the genetic backgrounds of the focal mouse and its Epha4 cage mate, as well as by sex. Observations of fighting in cages of males with dissimilar genetic backgrounds suggested aggression might underpin these modulatory effects, which prompted to us to expand the conceptual framework describing the mechanisms of social epistasis. Thus, our results strengthen the emerging view that IGE are not fixed properties of social genotypes but can depend on the genetic and sex context of interacting individuals, revealing an underexplored layer of their genetic architecture.