The evolution of religiosity by kin selection

ABSTRACT

Despite religion’s apparent ubiquity, hypotheses about the selection pressures that may have shaped its cognitive foundations remain controversial. Here, we develop and analyze a mathematical model inspired by Crespi and Summers’ suggestion that parent-offspring conflict has driven the evolution of religious beliefs to explore the causes and consequences of these selection pressures. To this end, we employ kin selection methodology to investigate how selection may mold an individual’s propensity for religiosity and corresponding patterns of gene expression, revealing that the evolution of religiosity is modulated by genetic relatedness between social partners, that selection in relation to religiosity may depend on an individual’s age and sex, and that religiosity can foment intragenomic conflicts of interest that give rise to parent-of-origin specific patterns of gene expression and concomitant clinical disorders. More generally, we develop a formal, theoretical framework that enables the derivation of clear-cut, comparative predictions about adaptive as well as maladaptive religiosity phenotypes.

KEYWORDS:

• Evolution
• genomic imprinting
• inclusive fitness
• intragenomic conflict
• imprinting disorders
• kin selection
• neurodevelopmental disorders
• religion

Introduction

Religion is regarded as a “human universal”, meaning that it has been found in all known contemporary and historical human societies (Brown, 2000). It appears to have existed from at least the Upper Palaeolithic (Sterelny, 2018) and has evolved into a complex and culturally diverse phenomenon (Wilson, 2002). This variety has led to a lack of consensus on how to even define religion, but a useful working definition is that religion represents a shared system of beliefs and actions concerning supernatural agency (Barrett, 2000, p. 29). It is widely acknowledged that humans are predisposed to imagine and accept supernatural concepts as products of our normal cognition but—given its counterintuitive contents, associated costly commitments and other seemingly maladaptive effects—understanding why religion should feature so prominently in human populations remains a substantial problem for evolutionary biology and allied disciplines (Powell & Clarke, 2012; Sosis, 2009).

There are two major approaches to explaining religion from an evolutionary standpoint (Table 1). Proponents of “by-product” hypotheses suggest that religion arises as an incidental consequence of cognitive processes that have themselves been evolutionarily favored for reasons having nothing to do with religion (Atran, 2005; Atran & Henrich, 2010; Barrett, 2000; Barrett & Lanman, 2008; Boyer, 2003; Guthrie, 1980; Hinde, 1999). A potential candidate for one of these cognitive processes underlying religiosity is “Theory of Mind” (e.g., Crespi, 2016; Gervais, 2013; McKay & Whitehouse, 2014), which enables humans to infer the mental states of others and therefore serves to facilitate interactions in one’s social environment. This capacity usefully applies even to social partners who are not physically present, and by extension—without requiring the invocation of further, special adaptations—it is easy to see how humans might also apply this capacity to deceased individuals and even non-animate entities, forming a basis for religiosity (see McKay & Whitehouse, 2014). However, the role of mentalizing in the expression of religiosity remains unclear (Di Dio et al., 2018; Ishii & Watanabe, 2021; Jack et al., 2016; Kapogiannis et al., 2009; Lindeman et al., 2015; Maji et al., 2017; Norenzayan et al., 2012; Reddish et al., 2016; Willard & Norenzayan, 2013).

Table 1. Religion as adaptation versus by-product. A classification of hypotheses on the evolution of religion according to mode of inheritance, unit of adaptation, and adaptive function—if any.

By-Product/Adaptation	Inheritance mode	Unit of adaptation	Adaptive function	References
By-Product	Genetic	Individual	None beyond primary functions of ordinary cognitive features	Guthrie, 1980; Hinde, 1999; Barrett, 2000; Boyer, 2001; Boyer, 2003; Atran, 2005; Dawkins, 2006; Dennet, 2006; Barrett & Lanman, 2008; Atran & Henrich, 2010
Adaptation	Genetic	Individual	Fostering intra-family cooperation	Crespi & Summers, 2014; Crespi, 2016
			Fostering (large-scale) intra-group cooperation	Wilson, 2002; Sosis & Alcorta, 2003; Alcorta & Sosis, 2005; Purzycki & Sosis, 2009; Bulbulia, 2008; Sosis, 2009; Bulbulia & Frean, 2010; Purzycki & Sosis, 2013; Powell & Clarke, 2012; Kiper & Sosis, 2014; Norenzayan et al., 2016
	Cultural	Group / Institution	Fostering (large-scale) intra-group cooperation	Wilson, 2002; Bulbulia, 2008; Atran & Henrich, 2010; Bulbulia & Frean, 2010; Powell & Clarke, 2012; Crespi & Summers, 2014; Kiper & Sosis, 2014; Crespi, 2016; Norenzayan et al., 2016; Szocik, 2017
		Meme	Dissemination via appealing to human cognition	Boyer, 2001; Dawkins, 2006; Dennet, 2006

In contrast, “adaptationist” hypotheses posit the existence of a “religious phenotype” that has emerged as an adaptation in its own right (Purzycki & Sosis, 2013). Adaptationist hypotheses subdivide into those which view religion as a product of genetical evolution, and representing an adaptation on the part of the individual for the purpose of fostering cooperation within families and/or within wider society (Alcorta & Sosis, 2005; Bulbulia, 2008; Bulbulia & Frean, 2010; Crespi, 2016; Crespi & Summers, 2014; Purzycki & Sosis, 2009; Purzycki & Sosis, 2013; Sosis, 2009; Sosis & Alcorta, 2003), versus those which view religion as a product of cultural evolution, and representing an adaptation of the social group/cultural institution (Atran & Henrich, 2010; Bulbulia, 2008; Bulbulia & Frean, 2010; Crespi, 2016; Crespi & Summers, 2014; Kiper & Sosis, 2014; Norenzayan et al., 2016; Powell & Clarke, 2012; Szocik, 2017; Wilson, 2002) or else the meme (Boyer, 2001; Dawkins, 2006; Dennet, 2006). The multifaceted nature of religion means that these hypotheses are not mutually exclusive with, for example, some seeing the biological capacity for religion as a simple by-product which has subsequently been hijacked by a cultural process driven by the evolutionary interests of memes (Dawkins, 2006).

Recently, it has been suggested that religion may have originated as a means to suppress intra-family conflict via parental manipulation. Crespi and Summers (Crespi, 2016; Crespi & Summers, 2014)—building upon and synthesizing the work of Alexander (2006), Lahti (2009), Coe and Palmer (2008, 2013), Palmer et al. (2008), Steadman and Palmer (2008), amongst others—have proposed that during human evolutionary history a mother may have been able to increase her inclusive fitness by instilling her children with beliefs in moralizing supernatural agents (e.g., deceased ancestors), thereby encouraging them to increase their cooperativeness with her and each other. They argue that children, owing to their strong dependence upon social learning, a tendency to readily adopt supernatural explanations and a pre-existing moral propensity, are predisposed to accept such manipulation. They have further suggested that religion may have subsequently been elaborated—both genetically and culturally—within and beyond the intra-family context as a means for social control and the establishment of cooperative relationships more generally. Following from this, they have predicted a tight linkage between the proximate mechanisms of religious and social bonding traits, placing a particular emphasis on the role of the neuropeptide oxytocin, which is understood to play a key part in the formation and maintenance of human social relationships, and which has been implicated in behaviors and cognitive features identified as important aspects of religion (Crespi, 2016; Crespi & Summers, 2014; and references therein). Crespi and Summers have also highlighted the possibility for developing comparative predictions linking inter-individual and inter-group variation in religious cognition and behaviors with variation in the strength and nature of within-group versus between-group competition, and variation in the benefits and costs of religious (i.e., mentalistic) versus mechanistic cognition.

This is a compelling origin story for religion. However, Crespi and Summers’ hypothesis remains underdeveloped. First, their argument has been developed in purely verbal terms, rather than formally in mathematical terms. Second, their parent-offspring conflict scenario applies very broadly, implying that maternal manipulation would be straightforwardly favored across a wide range of demographic and ecological contexts, such that it is difficult to see why quantitative variation in demographic and ecological factors would translate into variation in the intensity of religious beliefs and behavior, which limits the extent to which clear-cut, comparative predictions can be derived from their hypothesis. Third, although it is plausible that children would initially be susceptible to their mothers’ religious indoctrination it is unclear that this susceptibility would be evolutionarily maintained in the face of potentially strong selection to reduce their receptiveness for supernatural ideas, so as to avoid being manipulated.

Here, we undertake a formal treatment of the evolution of religion, building upon the verbal framework of Crespi and Summers. In contrast to Crespi and Summers, we focus our attention on the manipulated party—in the first instance, a child being exposed to supernatural ideas by their mother, and in the second instance, an adult being exposed to similar religious manipulation by their social partners—in order to investigate how natural selection shapes the individual’s susceptibility to such indoctrination. This involves describing the three-way tension between the direct-fitness Theory-of-Mind (and/or other) benefits associated with a cognition that predisposes the individual to supernatural ideas, the direct-fitness costs of being thereby more prone to manipulation, and the indirect-fitness (i.e., kin-selected) benefits arising from allowing oneself to be manipulated into being more cooperative with genetically-related social partners. The analysis enables the derivation of a suite of novel, concrete, comparative predictions concerning differences in religiosity in relation to age and sex, and also with respect to different elements of the genome, all modulated by variation in demographic and ecological factors. Consideration of intragenomic conflict with respect to religiosity further yields novel predictions concerning parent-of-origin specific patterns of expression in relation to genetic loci that affect religiosity and the maladaptive and clinical consequences of concomitant disorders of genomic imprinting. More generally, we intend our mathematical treatment to help connect a large and exciting literature on the possible evolutionary drivers of religion to current concepts and methodologies of social evolution theory.

Methods

We consider a large population divided into social groups, each containing women and men who are producing and raising children, with their success in this activity being modulated by religiosity—that is, an individual’s susceptibility to supernatural concepts. On the one hand, we assume that individuals with higher religiosity enjoy basic benefits of improved Theory-of-Mind and/or other abilities, which allow them to navigate through their social lives more successfully. On the other hand, we assume that individuals with higher religiosity tend to be less competitive with their social partners for resources—insofar as supernatural agents are invoked to manipulate one to behave more selflessly with social partners—which results in both a direct-fitness cost for themselves and an indirect-fitness benefit for their social partners, in line with the classic “tragedy of the commons” model (Frank, 1998; Hardin, 1968). For children, we assume that individuals compete with their maternal siblings for resources that improve their survival to adulthood. For adults, we assume that individuals compete with same-sex group mates for resources that improve their fecundity and investment into their children. That is, religiosity is associated with inclusive-fitness costs and benefits which we expect to balance out at equilibrium. For simplicity, we assume non-overlapping generations, with groups undergoing fission at the end of each generation and with daughter groups moving off to compete with other, unrelated groups for reproductive resources (Gardner & West, 2006; Haldane, 1932). Following density-dependent regulation that maintains the total number of groups in the population at a fixed level, adults disperse with sex-specific rates to other groups or else remain in their native group. Mating then occurs at random within each group, with some adults potentially achieving more reproductive success than their same-sex group mates, according to sex-specific degrees of reproductive skew. We perform a kin-selection analysis (Frank, 1998; Hamilton, 1964; Taylor, 1996; Taylor & Frank, 1996) to investigate how natural selection may act upon religiosity as a function of the individual’s age and sex, and under different assumptions of the genetic architecture of the trait (see Appendix for full details).

Results

Kin selection favors religiosity

Our model involves a three-way tension that modulates the evolution of religiosity. First, there is a direct-fitness benefit associated with greater predisposition to religious ideas owing to the cognitive processes that give rise to religiosity more generally, underpinning Theory of Mind and/or other abilities that are important for an individual’s success in navigating their social life. Second, a higher degree of religiosity is also associated with a greater susceptibility to manipulation by social partners, which makes the focal individual less competitive for resources. Third, this loss of competitiveness leads to an indirect-fitness (i.e., kin-selected) benefit owing to the increased success of the focal individual’s competitors, to the extent that they are their genetic relatives. Accordingly, in the absence of genetic relatedness the individual is favored to exhibit a degree of religiosity that exactly balances the direct-fitness Theory of Mind and/or other benefits against the direct-fitness costs due to loss of competitiveness, i.e., that which maximizes their direct fitness overall. And as social partners increasingly share genes in common with each other (higher relatedness), the individual is favored to exhibit a higher degree of religiosity than this baseline, resulting in an overall reduction in their direct fitness and an overall increase in the fitness of their social partners (Figure 1; Table 2, Prediction 1).

Figure 1. Kin selection favors religiosity. (a) In children, a higher degree of relatedness among maternal siblings (owing to greater probability of having the same father, β) favors a greater level of religiosity. (b) In adults, women (orange) favor a higher degree of religiosity due to higher within-group relatedness for women in scenarios with male-biased dispersal (d_M > d_F), and men (purple) favor a higher degree of religiosity due to higher within-group relatedness for men in scenarios with female-biased dispersal (d_F > d_M). (c) Women (orange) favor a higher degree of religiosity in scenarios with higher female reproductive variance (α > β), and men (purple) favor a higher degree of religiosity in scenarios with higher male reproductive variance (β > α).

Table 2. Predictions emerging from our analysis. All predictions are based on relatedness considerations and assume no sex differences in fitness costs and benefits of religiosity.

#	Prediction
1	Higher relatedness promotes religiosity.
2	In children, full (as opposed to half) sibship promotes religiosity, equally for both sexes.
3	Lower dispersal promotes religiosity in adulthood.
4	Higher sibship within social groups promotes religiosity in adulthood.
5	The least-dispersing sex has a higher potential for religiosity.
6	The sex with greater reproductive variance has a higher potential for religiosity.
7	Genomic imprinting will be favored at loci affecting religiosity.
8	In children, religiosity-promoter loci will be maternally expressed and paternally silenced, whereas religiosity-inhibitor loci will be paternally silenced and maternally expressed.
9	All else being equal, in adults, religiosity-promoter loci will be: paternally expressed and maternally silenced if dispersal is female-biased, and the reverse if dispersal is male-biased; paternally expressed and maternally silenced if reproductive variance is male-biased, and the reverse if reproductive variance is female-biased; and the opposite patterns will be true of religiosity-inhibiter loci.
10	Mutations and epimutations at religiosity loci will be associated with parent-of-origin specific maladaptive phenotypes in predictable directions, as detailed in Figure 3.

In line with Crespi and Summers (2014), we assume that religiosity in childhood modulates the individual’s cooperativeness towards—and, in turn, their competitiveness against—maternal siblings for maternal resources. Here, we consider that maternal siblings may or may not share the same father, and by varying the degree of their paternal sibship we are able to vary the degree of relatedness between them in order to investigate its impact on the evolution of religiosity in childhood. In accordance with the general prediction given above, we find that higher relatedness—owing to higher paternal sibship—leads to a greater level of religiosity being favored (Figure 1a). Note that, owing to the symmetrical inheritance of autosomal genes in relation to the two sexes, we do not expect there to be sex differences in the relatedness of maternal siblings. Hence, in the absence of sex differences in the costs and benefits of religiosity we expect there to be no sex differences in religiosity during childhood, with girls and boys being equally receptive to maternal religious manipulation (Table 2, Prediction 2).

In adults, we assume that religiosity modulates the individual’s competitiveness for reproductive resources against same-sex groupmates extending beyond the nuclear family. Accordingly, demographic and ecological factors that influence the degree of relatedness between groupmates are expected to impact upon the evolution of religiosity. We find that lower rates of dispersal and higher maternal and paternal sibship between group mates leads to a higher degree of within-group relatedness and hence a greater level of religiosity (Figure 1b & c; Table 2, Predictions 3 & 4).

Moreover, we find that sex differences in dispersal rate and/or unequal degrees of maternal and paternal sibship lead to sex differences in relatedness to group mates and hence to different levels of religiosity being favored in women and men, even if the basic costs and benefits of religiosity are exactly the same for both sexes. Specifically, because adults of the less-dispersing sex are, on average, more related to their group mates, then—all else being equal—in populations characterized by male-biased dispersal (e.g., matrilocality) religiosity is expected to be higher among women than among men, whereas in populations characterized by female-biased dispersal (e.g., patrilocality) religiosity is expected to be higher among men than among women (Figure 1b; Table 2, Prediction 5). Also, if reproduction by one of the sexes is dominated by a smaller number of individuals, then they will be more inclined towards religiosity because any group benefit will largely accrue to their own reproductive success, and accordingly—all else being equal—we would expect populations characterized by a higher degree of maternal sibship (e.g., polyandry) to favor higher religiosity among women than among men, whereas populations characterized by a higher degree of paternal sibship (e.g., polygyny) would favor higher religiosity among men than among women (Figure 1c; Table 2, Prediction 6).

Kin selection drives intragenomic conflict over religiosity

Above, we have shown that the evolution of religiosity is modulated by the degree of genetic relatedness between social partners. However, different parts of an individual’s genome may differ in the extent to which they share genes in common with the individual’s social partners. Crucially, individuals carry two genes at every one of their autosomal loci, one derived from the individual’s mother and one derived from the individual’s father, and these two genes are liable to be differently related to the individual’s social partners owing to sex differences in demographic and ecological factors (Haig, 1997). This means that genes deriving from the individual’s mother are liable to experience selection pressures in relation to the individual’s religiosity phenotype that are different from those experienced by the genes deriving from the individual’s father, and hence there is expected to be an intragenomic conflict (Burt & Trivers, 2006; Gardner & Úbeda, 2017; Haig, 1997) with each set of genes having a different optimum with regard to the individual’s level of religiosity (Figure 2).

Figure 2. Kin selection drives intragenomic conflict over religiosity. (a) In children (green), maternal-origin genes (dashed) always favor a higher degree of religiosity than paternal-origin genes (dotted) and genes that are ignorant of their origin (solid), except for in the special case whereby maternal offspring are guaranteed to have the same father. (b) In women (orange) and men (purple), maternal-origin genes (dashed) favor a higher degree of religiosity than paternal-origin genes (dotted) and genes that are ignorant of their origin (solid) in scenarios with male-biased dispersal (d_M > d_F), and the opposite pattern arises in scenarios with female-biased dispersal (d_F > d_M). (c) Maternal-origin genes (dashed) favor a higher degree of religiosity than paternal-origin genes (dotted) and genes that are ignorant of their origin (solid) in scenarios with higher female reproductive variance (α > β) and the opposite pattern arises in scenarios with higher male reproductive variance (β > α).

Such differences in gene interests are understood to drive the evolution of parent-of-origin specific gene expression, or “genomic imprinting” (Haig, 1997). Specifically, if a locus encodes a gene product that increases the individual’s religiosity (a “religiosity promoter”), then the gene with the higher optimum is expected to favor a greater degree of gene expression and its homologue is expected to favor a lower degree of gene expression, culminating in the silencing of the latter gene and the expression of the former gene at its optimal level, in what has been termed the “loudest voice prevails” principle (Haig, 1996). Conversely, if the gene product decreases the individual’s religiosity (a “religiosity inhibitor”) then the reverse pattern of gene expression is predicted. Accordingly, if kin selection has been a driver of the evolution of religiosity, we would expect religiosity loci to show a greater tendency towards parent-of-origin specific gene expression (Table 2, Prediction 7).

During childhood, the salient social partners are the individual’s maternal siblings, who may or may not be paternal siblings. This means that, except for in populations where all maternal siblings are guaranteed to share the same father, relatedness will generally be higher with respect to maternal-origin genes than with respect to paternal-origin genes, and hence we expect a child’s maternal-origin genes will favor a greater level of religiosity than will the paternal-origin genes (Figure 2a). Accordingly, during childhood, religiosity-promoter loci are expected to be maternally-expressed and paternally-silenced, whereas religiosity-inhibitor loci are expected to be paternally-expressed and maternally-silenced (Figure 3; Table 2, Prediction 8).

Figure 3. Genomic imprinting and maladaptation in relation to religiosity. Genomic imprinting is predicted to arise as a result of intragenomic conflict over religiosity, and to be associated with maladapted phenotypes including in relation to four classes of mutational perturbation (gene deletions, gene duplications, epimutations, and uniparental disomies). The patterns predicted for adult religiosity phenotypes in this table are for scenarios in which relatedness is higher for the paternal-origin gene rather than for the maternal-origin gene (i.e. opposite from the situation for children interacting socially with maternal siblings).

In adults, we also find potential for intragenomic conflict over religiosity. As in children, maternal-origin genes and paternal-origin genes in adults can favor different optimal levels of religiosity. During adulthood, however, relatedness between group members may depend upon dispersal as well as mating patterns and can therefore vary significantly across different populations. All else being equal, in populations with male-biased dispersal, relatedness will generally be higher with respect to maternal-origin genes than with respect to paternal-origin genes, and as a consequence we again expect an individual’s maternal-origin genes to favor a greater level of religiosity than will their paternal-origin genes (see left half of Figure 2b), such that religiosity-promoter loci are expected to be maternally-expressed and paternally-silenced and religiosity-inhibitor loci are expected to be maternally-silenced and paternally-expressed. Conversely, all else being equal, in populations with female-biased dispersal, we expect paternal-origin genes will favor a greater level of religiosity than will maternal-origin genes (see right half of Figure 2b), such that religiosity-promoter loci are expected to be maternally-silenced and paternally-expressed and religiosity-inhibitor loci are expected to be maternally-expressed and paternally-silenced (Figure 3).

All else being equal, in polyandrous populations (with higher maternal than paternal sibship) we expect an adult’s maternal-origin genes to favor a greater level of religiosity than will their paternal-origin genes (see left half of Figure 2c), such that religiosity-promoter loci are expected to be maternally-expressed and paternally-silenced and religiosity-inhibitor loci are expected to be maternally-silenced and paternally-expressed, and in polygynous populations (with higher paternal than maternal sibship) we expect an adult’s paternal-origin genes to favor a greater level of religiosity than their maternal-origin genes (see right half of Figure 2c), such that religiosity-promoter loci are expected to be maternally-silenced and paternally-expressed and religiosity-inhibitor loci are expected to be maternally-expressed and paternally-silenced (Figure 3). Interestingly, in those scenarios involving within-group relatedness being higher via paternal-origin genes, we expect the pattern of imprinting of religiosity loci to reverse between childhood and adulthood, from maternally-expressed and paternally-silenced to maternally-silenced and paternally-expressed at promoter loci and the other way round at inhibitor loci (Figure 3; Table 2, Prediction 9).

Since only one of the two gene copies at an imprinted locus is expressed, the individual is functionally haploid at this locus and therefore potentially more vulnerable to the detrimental effects of mutation (Wilkins & Haig, 2003). For example, a loss-of-function mutation in the expressed gene would result in the complete absence of the functional gene product at the imprinted locus, with possibly drastic phenotypic consequences, as opposed to a mere halving of the functional gene product that would be expected at a non-imprinted locus. Moreover, imprinting may result in mutations having different phenotypic effects according to their parent of origin. For example, the above loss-of-function mutation would be expected to have no impact upon the phenotype if it occurred in the non-expressed gene at the same locus. We can therefore generate additional predictions about patterns concerning maladaptive phenotypes resulting from the evolution of religiosity (Figure 3). For instance, if we consider a promoter locus for religiosity in children, we expect the paternal-origin gene to be silenced and the maternal-origin gene to be expressed at a level associated with a “normal” religiosity phenotype. Accordingly, mutational deletion of the maternal-origin gene is expected to result in the complete absence of the religiosity-promoting gene product, and therefore a “hyporeligious” phenotype. In contrast, deletion of the paternal-origin gene is expected to have no effect, i.e., giving rise to a “normal” phenotype. The parent-of-origin specific phenotypic consequences of a suite of other mutational and epimutational perturbations can be similarly determined (Figure 3; Table 2, Prediction 10).

Discussion

Despite the apparent ubiquity of religion across human populations, the selection pressures that have shaped its cognitive foundations have remained obscure. Crespi and Summers (2014; see also Crespi, 2016) have suggested that the origins of religion might lie in parent-offspring conflict, whereby mothers employed religious indoctrination as means of manipulating their children into more cooperative behavior, and that religion has subsequently spread to become a manipulative tool employed more generally in the context of social conflict. We have developed and analyzed a mathematical, kin-selection model inspired by this hypothesis. By refocusing attention on the manipulated party, we have derived a suite of new comparative predictions concerning variation in religiosity as a function of sex, age and ecological context, and also parent-of-origin specific patterns of gene expression for loci underpinning religiosity and concomitant maladaptive phenotypes associated with a range of mutational perturbations, which present novel avenues for empirical testing across multiple disciplinary domains.

Our analysis has revealed that kin selection can promote the evolution of religiosity, with an individual’s predisposition to religious ideas being an increasing function of relatedness to social partners, on account of the benefits to kin offsetting some of the personal costs of religious manipulation. This strengthens the plausibility of Crespi and Summers’ (2014) hypothesis, by helping to explain why individuals would not be favored to reduce their vulnerability to indoctrination, and also yields new comparative predictions concerning how variation in religiosity across populations is characterized by different demographies. In relation to maternal manipulation, we have found that children should be more readily accepting of religious indoctrination that promotes cooperation with their maternal siblings as the likelihood of their sharing the same father increases, on account of this being associated with greater relatedness between siblings. Following from this, we would expect that children in traditionally more monogamous and polygynous populations to be more predisposed to religiosity than children in polyandrous populations. Similarly, we have found that the susceptibility of adults to religious indoctrination should also be higher when they are more related to their social partners, as for example in populations characterized by lower rates of dispersal and higher degrees of reproductive skew.

Our analysis also yields new comparative predictions concerning how religiosity varies within populations—in particular, between the sexes. In the context of our model there are no relatedness differences between girls and boys with respect to their maternal siblings, which gives no basis for expecting sex differences in religiosity to manifest in childhood. However, if the sexes disperse at different rates then we expect them to experience different degrees of relatedness to social partners in adulthood, with individuals of the least-dispersing sex tending to interact with more highly related social partners and hence being favored to have a greater predisposition to religiosity. For example, if men predominantly move to other groups in order to live with their spouse and the latter’s kin (matrilocality), as has been suggested for ancestral Austronesian societies (Jordan et al., 2009), we would expect women, who tend to stay in their natal groups, to be more related to each other than are men, and hence that religiosity would be higher in women. In a population with a polyandrous mating system, we would expect that this effect would be even stronger. Phylogenetic analyses from contemporary hunter-gatherer societies, however, indicate that serial monogamy or low-level polygyny are more likely mating systems for ancestral populations (Walker et al., 2011). If this was the case for ancestral Austronesian societies, reproductive variance would be higher for men than women with the consequence that we would expect that men would favor a higher degree of religiosity compared to women.

In addition, we have found potential for within-individual, intragenomic conflict between maternal-origin and paternal-origin genes in both children and adults, to the extent that there is parent-of-origin information available. Differential within-group relatedness via maternal-origin versus paternal-origin genes leads to these genes favoring different levels of religiosity and may consequently lead to the evolution of parent-of-origin specific gene expression, i.e., “genomic imprinting” (Haig, 1997). By considering loci that either promote or inhibit the expression of traits associated with religiosity, our predictions regarding imprinting status for different phenotypes (see Figure 3) can be tested directly and furthermore, may give us insight into historical dispersal and mating patterns. In populations in which children are raised alongside maternal siblings who need not be paternal siblings, their maternal-origin genes are expected to favor a higher degree of religiosity than their paternal-origin genes, such that religiosity-promoter loci are expected to be maternally-expressed and paternally-silenced and religiosity-inhibitor loci are expected to be maternally-silenced and paternally-expressed, except in the complete absence of female promiscuity (i.e., maternal siblings always have the same father). In contrast, we expect patterns of genomic imprinting in adults to be more dependent on sex-specific patterns of dispersal and reproductive skew. Indeed, in populations characterized by female-biased dispersal (i.e., patrilocality) and/or male-biased reproductive skew (i.e., polygyny) we expect an adult’s paternal-origin genes to favor a higher level of religiosity than their maternal-origin genes, and hence a pattern of genomic imprinting exactly opposite to that for children in the same population. This implies that, for any locus whose action influences religiosity both in children and in adults, the pattern of imprint will reverse as the individual ages from childhood to adulthood.

Our analysis also yields new predictions concerning maladaptive phenotypes arising from a range of mutational perturbations. These predictions provide an avenue for improved understanding of clinical disorders manifesting a religiosity dimension. For instance, mis-expression of imprinted genes can lead to phenotypically diametric disorders such as Beckwith-Wiedemann and Silver-Russell Syndromes, or Prader-Willi and Angelman Syndromes, with profound effects on pre- and post-natal growth, adult metabolism, and social cognition (Ishida & Moore, 2013; Kalish et al., 2014; Millership et al., 2019; Peters, 2014; Plasschaert & Bartolomei, 2014; Wilkinson et al., 2007). To the extent that religiosity is associated with mentalistic cognition, it is worth considering disorders that affect Theory-of-Mind and related capabilities, such as autism and psychosis (e.g., Gray et al., 2011; also see Hill & Frith, 2003): whilst there is evidence of a positive association between schizotypal traits and aspects of religiosity (e.g., Barnes & Gibson, 2013; Breslin & Lewis, 2015; Iyassu et al., 2014; Lindeman & Lipsanen, 2016), the relationship between autistic traits and religiosity appears more complex (Jack et al., 2016; Lindeman & Lipsanen, 2016; Norenzayan et al., 2012; Reddish et al., 2016). Crespi and Badcock (2008) have argued that autism and psychosis represent opposite extremes of mentalistic-mechanistic cognition (see e.g., Thakkar et al., 2008, for an opposing view), and that this has a bearing on the clinical consequences of genomic-imprinting disorders. This “imprinted-brain” theory holds that negative symptoms in autism-spectrum disorders, and positive symptoms in psychotic-spectrum conditions, represent hypo-mentalistic, or paternally biased, and hyper-mentalistic, or maternally biased cognition, respectively, resulting from dysregulated imprinting. For example: lack of expression of the maternally expressed gene UBE3A in the 15q11-q13 chromosomal region is implicated in the pathogenesis of Angelman Syndrome, which exhibits autism-spectrum characteristics; overexpression of UBE3A is implicated in the pathogenesis of Prader-Willi Syndrome, which exhibits psychosis-related characteristics; and genetic variation in UBE3A architecture is associated with variation in total schizotypy, i.e., degree of psychosis-related characteristics (Salminen et al., 2019). Since religious delusions are prevalent in psychosis-spectrum disorders (Anderson-Schmidt et al., 2019), the 15q11-q13 chromosomal region presents a possible focus for future investigation into genetic influences upon religiosity and—combined with data on ancestral mating and dispersal patterns—empirical testing of our predictions.

Our aim has been to investigate how patterns of genetic relatedness translate into clear-cut, comparative predictions that may serve to illuminate the selective forces that have modulated the evolution of religiosity. To this end, we have not explicitly considered the consequences of individual variation in the costs and benefits. This shortcoming is particularly acute in relation to our predictions concerning sex differences in religiosity, as women and men are liable to experience different personal-fitness consequences of their own and their social partners’ religiosity phenotypes, and such effects are liable to confound the patterns we have described here. A general driver of sex differences in costs and benefits is sexual selection, which across many species—including humans—is expected to operate differently between females and males (Andersson, 1994; Darwin, 1871), and indeed sexual selection has been suggested to be a key driver of the evolution of religion itself (Miller, 2007; see also Soler & Lenfesty, 2016; for an in depth summary of hypotheses see Czachesz, 2018). An exploration of sex-specific costs and benefits represents a major avenue for future exploration on this topic. In contrast, our predictions concerning intragenomic conflicts and concomitant patterns of imprinting and associated clinical disorders are expected to be robust to sex differences in costs and benefits as they hinge upon inclusive-fitness differences between genes that reside in the very same bodies, and thereby naturally control for such confounding variables (cf. Rautiala & Gardner, 2016).

Similarly, although the scenario described by Crespi and Summers (2014) and investigated here involves a flow of ideas concerning supernatural agents from mother to child, and between group mates more generally, our explicit focus has been on the genetical molding of the religiosity phenotype and accordingly the cultural dynamics of these religious beliefs have remained implicit. Our analysis seeks to understand the consequences of there being at least some genetic variation in religiosity, and it does not require that all—or even most—of the variation in such traits has a genetic basis. In combination with research into ancestral ecologies, our predictions as to the relationships between religiosity and ecological parameters yield insights not only into the selective origins of religion, but also into the demographic circumstances that influenced its expression. Nonetheless, we anticipate cultural transmission dynamics to also play a crucial role in determining these traits and systems, including in relation to gene-culture co-evolution (e.g., Bulbulia, 2008; Czachesz, 2018; Ferretti & Adornetti, 2014; Norenzayan et al., 2016; Rowthorn, 2011; Szocik, 2017). For example, it has been suggested that ancestor worship has evolved as a descendant-leaving strategy with culturally learned, cross-generationally transmitted cooperation among descendants—incentivized by supernatural concepts—leading to greater success for one’s lineage (Clark & Coe, 2021; Coe et al., 2010). Such an approach could represent a bridge between our kin-selection approach and more culturally focused approaches to investigating the evolution of religious systems. Moreover, the possible role for religion to modulate how individuals from different backgrounds with different religious systems interact with each other might also have had a major impact upon population demography and genetic evolution. These represent exciting avenues for future exploration.

Acknowledgements

We thank Carolina Barata, Mauricio González-Forero, Thomas Hitchcock, Jasmeen Kanwal, Petri Rautiala and an anonymous reviewer for useful comments and discussion.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was supported by a Natural Environment Research Council Independent Research Fellowship (grant number NE/K009524/1) and a H2020 European Research Council Consolidator (grant number 771387).

Appendix

General analysis of religiosity

We can express a focal juvenile’s relative fitness W in terms of their own investment into religiosity (x_FJ if they are female and x_MJ if they are male), their siblings’ average investment into religiosity (y_FJ for their sisters and y_MJ for brothers), their parents’ investment into religiosity (x_FA for their mother and x_MA for their father), their parents’ social partners’ average investment into religiosity (y_FA for women and y_MA for men), as well as the population averages of these quantities (z_FJ, z_MJ, and z_FA, z_MA), where F and M denote female and male, and J and A denote juvenile and adult, respectively. We assume that the fitness function is symmetrical with respect to the sex of actor and recipient, such that the only possible sex differences are with respect to the phenotypes that individuals express.

We consider an autosomal locus G_ij which influences the degree of religiosity of a class-ij individual, where i ∈ {F, M} and j ∈ {J, A}. Drawing one of the focal individual’s two genes at random from this locus and denoting its genic value by g_ij, then natural selection acts to increase the population average genic value—and hence the average investment made by class-ij individuals into religiosity—if dW/dg_ij > 0, where the derivative is evaluated at the population average (Taylor & Frank, 1996). Since a given carrier may be female or male, relative fitness is given as a weighted average taken across female and male juveniles, i.e., W = c_F W_F + c_M W_M (Taylor, 1996), where c_F = c_M = ½ are the class reproductive values of female and male juveniles, respectively (Fisher, 1930; Price, 1970; Taylor, 1996). Assuming that a gene’s impact on the phenotype does not depend on its parent of origin, we can rewrite the left-hand side of the above condition, using the chain rule, as (1) ${\displaystyle \frac{\mathrm{d}W}{\mathrm{d}{g}_{\mathrm{i}\mathrm{j}}}=\sum _{\mathrm{k}\in \left\{\mathrm{F},\mathrm{M}\right\}}{c}_{\mathrm{k}}\sum _{\begin{array}{c}\mathrm{l}\in \left\{\mathrm{F},\mathrm{M}\right\}\\ \mathrm{m}\in \left\{\mathrm{J},\mathrm{A}\right\}\end{array}}\left({\displaystyle \frac{\mathrm{\partial }{W}_{\mathrm{k}}}{\mathrm{\partial }{x}_{\mathrm{l}\mathrm{m}}}{\displaystyle \frac{\mathrm{\partial }{x}_{\mathrm{l}\mathrm{m}}}{\mathrm{\partial }{G}_{\mathrm{l}\mathrm{m}}}{\displaystyle \frac{\mathrm{d}{G}_{\mathrm{l}\mathrm{m}}}{\mathrm{d}{g}_{\mathrm{k}\mathrm{i}\mathrm{j}}}+{\displaystyle \frac{\mathrm{\partial }{W}_{\mathrm{k}}}{\mathrm{\partial }{y}_{\mathrm{l}\mathrm{m}}}{\displaystyle \frac{\mathrm{\partial }{y}_{\mathrm{l}\mathrm{m}}}{\mathrm{\partial }{G}_{\mathrm{l}\mathrm{m}}^{\mathrm{\prime }}}{\displaystyle \frac{\mathrm{d}{G}_{\mathrm{l}\mathrm{m}}^{\mathrm{\prime }}}{\mathrm{d}{g}_{\mathrm{k}\mathrm{i}\mathrm{j}}}}}}}}}\right),}$(1) where: G_lm is the focal individual's genetic value (if m = J) or the focal individual’s sex-l parent’s genetic value (if m = A), for the class-lm investment into religiosity; G_lm’ is the average genetic value of the focal individual's sex-l siblings (if m = J) or the focal individual’s parents’ sex-l social partners (if m = A), for the class-lm investment into religiosity; ∂x_lm/∂G_lm = ∂y_lm/∂G_lm’ = γ_lm represents the mapping of class-lm religiosity genetic value to class-lm religiosity phenotype; ∂G_lm/∂g_kij represents the genetic association between a sex-k juvenile’s class-ij religiosity genic value and either their own (if m = J) or their sex-l parent’s (if m = A) class-lm religiosity genetic value; and dG_lm’/dg_kij represents the genetic association between a sex-k juvenile’s class-ij religiosity genic value and either their sex-l siblings’ (if m = J) or their parents’ sex-l social partners’ (if m = A) class-lm religiosity genetic value, where l ∈ {F, M} and m ∈ {J, A}. Note that ∂W_k/∂x_lm = 0 if l ≠ k and m = J, as the focal juvenile’s fitness is not a function of the phenotype it would have expressed had it been a member of the opposite sex. Also note that if l = i and m = j then dG_lm/dg_kij = p_kij represents the consanguinity of a sex-k juvenile to themselves (if m = j = J) or to their sex-l parent (if m = j = A) and dG_lm’/dg_kij = p_kij’ represents the consanguinity of a sex-k juvenile to their sex-l siblings (if m = j = J) or to their parents’ sex-l social partners (if m = j = A), and if l ≠ i and/or m ≠ j then dG_lm/dg_kij = dG_lm’/dg_kij = 0 upon the assumption that there is no pleiotropy or linkage disequilibrium between the different classes’ religiosity traits. Further note that, owing to the symmetries of diploid inheritance, p_kij = p_ij and p_kij’ = p_ij’ for all k ∈ {F,M}. Accordingly, natural selection favors an increase in religiosity, if (2) $-C\left(z_{\mathrm{i}\mathrm{j}}\right)+B\left(z_{\mathrm{i}\mathrm{j}}\right)r_{\mathrm{i}\mathrm{j}}>0,$(2) where: r_ij = p_ij’/ p_ij is the kin-selection coefficient of relatedness between the actor and recipient, which is determined by the respective consanguinities which may differ according to dispersal patterns and mating system (see below; Bulmer, 1994); -C(z_ij) = c_F ∂W_F/∂x_lm + c_M ∂W_M/∂x_lm is the marginal direct fitness effect of increased investment into the trait by the focal individual during childhood or by the parents during adulthood respectively; and B(z_ij) = c_F ∂W_F/∂y_lm + c_M ∂W_M/∂y_lm is the marginal indirect fitness cost or benefit of increased investment into the trait by the juvenile maternal siblings or by the adult group members, respectively.

Note that p_ij and p_ij^’ are strictly speaking the consanguinities of the recipient to the part of the relevant actor’s genotype that controls the actor’s phenotype. If both genes at a locus share equal control over the actor’s phenotype, then we have p_ij = p_{ij |I} = 1/2 p_{ij |M} + 1/2 p_{ij |P} and p_ij’ = p_{ij |I}’ = 1/2 p_{ij |M}’ + 1/2 p_{ij |P}’, where p_{ij |I} represents the consanguinity of a juvenile to themselves (if j = J) or to their sex-I parent (if j = A), p_{ij |M} represents the juvenile’s consanguinity to this individual’s maternal-origin gene, p_{ij |P} represents the juvenile’s consanguinity to this individual’s paternal-origin gene, p_{ij |I}’ represents the consanguinity of the juvenile to their sex-i siblings (if j = J) or to their parents’ sex-i social partners (if j = A), p_{ij |M}’ represents the juvenile’s consanguinity to these individuals’ maternal-origin genes, p_{ij |P}’ represents the juvenile’s consanguinity to these individuals’ paternal-origin genes, and relatedness is given by r_ij = r_{ij |I} = p_{ij |I}’/p_{ij |I}. If instead the actor’s phenotype is fully controlled by only their maternal-origin gene, then we have p_ij = p_ij|M and p_ij’ = p_ij|M’, and relatedness is given by r_ij = r_{ij |M} = p_{ij |M}’/p_{ij |M}. And if the actor’s phenotype is fully controlled by only their paternal-origin gene, then we have p_ij = p_{ij |P} and p_ij’ = p_{ij |P}’, and relatedness is given by r_ij = r_{ij |P} = p_{ij |P}’/p_{ij |P}.

The condition for increase takes the form -C(z_ij) + B(z_ij) ρ > 0, where ρ = r_{ij |M} for maternal-origin control, ρ = r_{ij |P} for paternal-origin control, and ρ = r_ij for equal control. Assuming that an intermediate, convergence-stable equilibrium z_ij* exists, then we can define a function J(z_ij*, ρ) = -C(z_ij*) + B(z_ij*) ρ such that J(z_ij*, ρ) = 0 and ∂J/∂z_ij* < 0 (Christiansen, 1991; Taylor, 1996). Using the chain rule, we can write dJ/dρ = (∂J/∂ρ) + (∂J/∂z_ij*)(dz_ij*/ dρ) = 0, which rearranges as dz_ij*/dρ = -(∂J/∂ρ)/ (∂J/∂z_ij*), and hence S(dz_ij*/dρ) = S(∂J/∂ρ) = S(B(z_ij*)) where the function S returns the sign of its argument (positive, negative or zero) (Farrell et al., 2015; Pen, 2000). Consequently, if the religiosity of social partners improves the focal individual’s fitness (B > 0), then higher relatedness is associated with a higher religiosity optimum (dz_ij*/dρ > 0); if the religiosity of social partners decreases the focal individual’s fitness (B < 0), then higher relatedness is associated with a lower religiosity optimum (dz_ij*/dρ < 0); and if the religiosity of social partners does not affect the focal individual’s fitness (B = 0), then higher relatedness is not associated with a higher or lower religiosity optimum (dz_ij*/dρ = 0).

Therefore, assuming that the religiosity of social partners improves the focal individual’s fitness (B > 0), then: higher relatedness is associated with a higher religiosity optimum (dz_ij*/dρ > 0), leading to results 1–4 of the main text; the religiosity optimum is higher for women than it is for men (z_FA* > z_MA*) if relatedness is higher for women than men (r_FA > r_MA), and the religiosity optimum is lower for women than for men (z_FA* < z_MA*) if relatedness is lower for women than men (r_FA < r_MA), leading to results 5 and 6 of the main text; maternal-origin genes will favor a higher religiosity optimum than paternal-origin genes (z_{ij |M}* > z_{ij |I}* > z_{ij |P}*) when relatedness is higher for the former than the latter (r_{ij |M} > r_ij > r_{ij |P}) and maternal-origin genes will favor a lower religiosity optimum than paternal-origin genes (z_{ij |M}* < z_{ij |I}*< z_{ij |P}*) when relatedness is lower for the former than for the latter (r_{ij |M} < r_{ij |I} < r_{ij |P}), leading to results 7–9 of the main text. If the religiosity of social partners decreases the focal individual’s fitness (B < 0) then our predictions are exactly reversed, and if the religiosity of social partners does not affect the fitness of the focal individual (B = 0), then the religiosity optimum for females is equal to that for males (z_FA* = z_MA*), and the religiosity optimum for maternal-origin genes is equal to that for genes ignorant of their origin as well as paternal-origin genes (z_{ij |M}* = z_{ij |I}* = z_{ij |P}*), with none of these quantities being dependent upon relatedness.

Illustration

To illustrate how this general analysis applies to hypotheses about specific functions of religiosity, we construct a simple model which incorporates associated cognitive properties of religiosity in juveniles, adults, females, and males, in a variety of demographic scenarios. We assume an even sex ratio at birth, such that K girls and K boys are born in each group, where K is a large constant. We denote the probability that two juveniles born in the same group share the same mother by α and the probability that they share the same father by β, and we assume that every child is the product of an independent, random pairing of a woman and man such that the probability that two maternal siblings share the same father is also β; this allows us to explore the effects of modulating the degree of polyandry (by varying α) and the degree of polygyny (by varying β). We denote the probability of a child’s survival to adulthood by S_k, and we assume that this is a function of their own (x_kJ), their parents’ (x_FA and x_MA), their maternal siblings’ (y_FJ and y_MJ) and their parents’ group mates’ (y_FA, y_MA) investment into religiosity. All adults of the parental generation then die, such there is no overlapping of generations. The surviving individuals within each group organize themselves into smaller groups or “buds” at random with some of their peers—with all buds containing the same number of individuals, and having an even sex ratio, and with any excess of individuals of one sex that are not incorporated into buds being assumed to perish—with the buds then dispersing to random locations elsewhere in the population (i.e., group fissioning; Gardner & West, 2006) and competing with the other groups that have also dispersed there for control of the resources in that location. One group is chosen at random to survive this competition in each location to be the parents of the next generation of children to be born there. Following this density-dependent regulation that maintains the total number of groups in the population at a fixed level, adults disperse with sex-specific rates—d_F for women and d_M for men—to other groups, occupying places vacated by other dispersers, or else remain in their native group.

We assume that a juvenile’s survival is given by the product of (i) a theory-of-mind benefit equal to their own level of religiosity, (ii) their share of parental resources obtained in competition with their maternal siblings; (iii) their mother’s share of group reproductive resources obtained in competition with other women; and (iv) their father’s share of reproductive resources obtained in competition with other men. We assume that resources are shared according to a “tragedy of the commons” (Frank, 1998) scenario, whereby the proportion of resources seized by an individual is proportional to their competitiveness and the total amount of resources is proportional to the extent that the group members refrain from outright competitiveness, and that an individual’s competitiveness is equal to one unit minus their level of religiosity. Accordingly, a sex-k juvenile survives to adulthood with probability (3) $S_k=x_{\mathrm{k}\mathrm{J}}{\displaystyle \frac{1-x_{\mathrm{k}\mathrm{J}}}{1-\left({\displaystyle \frac{1}{2}{y}_{\mathrm{F}\mathrm{J}}+{\displaystyle \frac{1}{2}{y}_{\mathrm{M}\mathrm{J}}}}\right)}\left({\displaystyle \frac{1}{2}{y}_{\mathrm{F}\mathrm{J}}+{\displaystyle \frac{1}{2}{y}_{\mathrm{M}\mathrm{J}}}}\right)x_{\mathrm{F}\mathrm{A}}{\displaystyle \frac{1-x_{\mathrm{F}\mathrm{A}}}{1-y_{\mathrm{F}\mathrm{A}}}\left({\displaystyle \frac{1}{2}{y}_{\mathrm{F}\mathrm{A}}+{\displaystyle \frac{1}{2}{y}_{\mathrm{M}\mathrm{A}}}}\right)x_{\mathrm{M}\mathrm{A}}{\displaystyle \frac{1-x_{\mathrm{M}\mathrm{A}}}{1-y_{\mathrm{M}\mathrm{A}}}\left({\displaystyle \frac{1}{2}{y}_{\mathrm{F}\mathrm{A}}+{\displaystyle \frac{1}{2}{y}_{\mathrm{M}\mathrm{A}}}}\right),}}}$(3) A sex-k juvenile’s expected relative fitness is W_k = S_k/ S¯_k if sex-k is the rarer sex and is W_k = (S_k/S′_k) × (S′_­k/S¯_­k) if sex-k is the more-common sex, where S′_k is the average survival of sex-k juveniles in the focal individual’s group, S′_­k is the average survival of non-sex-k juveniles in the focal individual’s group, S¯_k is the average survival of sex-k juveniles across the whole population, and S¯_­k is the average survival of non-sex-k juveniles across the whole population. Also note that: for a juvenile trait (j = J) we have p_ij = q_self and p_ij’ = q_sib; for an adult trait (j = A) we have p_ij = q_par and p_ij’ = q_F (if i = F) or p_ij’ =  q_M (if i = M); and that r_sib = q_sib/q_self, r_F = q_F/q_par and r_M = q_M/q_par, where r_sib = r_{sib |I} = q_{sib |I} /q_{self |I}, r_F = r_{F |I} q_{F |I} /q_{par |I}, and r_M = r_{M |I} = q_{M |I} /q_{par |I}, if both genes at a locus share equal control over the actor’s phenotype, r_sib = r_{sib |M} = q_{sib |M} /q_{self |M}, r_F = r_{F |M} q_{F |M} /q_{par |M}, and r_M = r_{M |M} = q_{M |M} /q_{par |M}, if the actor’s phenotype is fully controlled by only their maternal-origin gene, and r_sib = r_{sib |P} = q_{sib |P} /q_{self |P}, r_F = r_{F |P} q_{F |P} /q_{par |P}, and r_M = r_{M |P} = q_{M |P} /q_{par |P}, if the actor’s phenotype is fully controlled by only their paternal-origin gene.

Evaluating the cost and benefit terms in expression (2) for this illustrative fitness function, we obtain the condition for natural selection to favor an increase in the level of religiosity exhibited by female juveniles as (4) ${\displaystyle \frac{1-2z_{\mathrm{F}\mathrm{J}}}{2z_{\mathrm{F}\mathrm{J}}\left(1-z_{\mathrm{F}\mathrm{J}}\right)}+{\displaystyle \frac{2}{\left(z_{\mathrm{F}\mathrm{J}}+z_{\mathrm{M}\mathrm{J}}\left)\right(2-z_{\mathrm{F}\mathrm{J}}-z_{\mathrm{M}\mathrm{J}}\right)}{r}_{\mathrm{s}\mathrm{i}\mathrm{b}}+{ϵ}_{\mathrm{F}\mathrm{J}}>0,}}$(4) where ϵ_FJ = [(1-2z_FJ)/(2z_FJ(1-z_FJ))]r_group if z_FJ(1-z_FJ) > z_MJ(1-z_MJ) and ϵ_FJ = -[(1-2z_FJ)/(2z_FJ(1-z_FJ))]r_group if z_FJ(1-z_FJ) < z_MJ(1-z_MJ), and where r_group is the relatedness of two juveniles born in the same group. Similarly, we find that the condition for natural selection to favor an increase in the level of religiosity exhibited by male juveniles is (5) ${\displaystyle \frac{1-2z_{\mathrm{M}\mathrm{J}}}{2z_{\mathrm{M}\mathrm{J}}\left(1-z_{\mathrm{M}\mathrm{J}}\right)}+{\displaystyle \frac{2}{\left(z_{\mathrm{F}\mathrm{J}}+z_{\mathrm{M}\mathrm{J}}\left)\right(2-z_{\mathrm{F}\mathrm{J}}-z_{\mathrm{M}\mathrm{J}}\right)}{r}_{\mathrm{s}\mathrm{i}\mathrm{b}}+{ϵ}_{\mathrm{M}\mathrm{J}}>0,}}$(5) where ϵ_MJ = -[(1-2z_MJ)/(2z_MJ(1-z_MJ))]r_group if z_FJ(1-z_FJ) > z_MJ(1-z_MJ) and ϵ_MJ = [(1-2z_MJ)/(2z_MJ(1-z_MJ))]r_group if z_FJ(1-z_FJ) < z_MJ(1-z_MJ), the condition for natural selection to favor an increase in the level of religiosity exhibited by women is (6) ${\displaystyle \frac{\left(z_{\mathrm{F}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}\left)\right(1-2z_{\mathrm{F}\mathrm{A}}\right)}{2z_{\mathrm{F}\mathrm{A}}\left(1-z_{\mathrm{F}\mathrm{A}}\right)\left(z_{\mathrm{F}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}\right)}+{\displaystyle \frac{z_{\mathrm{F}\mathrm{A}}\left(2-z_{\mathrm{F}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}\right)}{2z_{\mathrm{F}\mathrm{A}}\left(1-z_{\mathrm{F}\mathrm{A}}\right)\left(z_{\mathrm{F}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}\right)}{r}_{\mathrm{F}}>0,}}$(6) and the condition for natural selection to favor an increase in the level of religiosity exhibited by men is (7) ${\displaystyle \frac{\left(z_{\mathrm{F}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}\left)\right(1-2z_{\mathrm{M}\mathrm{A}}\right)}{2z_{\mathrm{M}\mathrm{A}}\left(1-z_{\mathrm{M}\mathrm{A}}\right)\left(z_{2\mathrm{F}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}\right)}+{\displaystyle \frac{z_{\mathrm{M}\mathrm{A}}\left(2-z_{\mathrm{M}\mathrm{A}}+z_{\mathrm{F}\mathrm{A}}\right)}{2z_{\mathrm{M}\mathrm{A}}\left(1-z_{\mathrm{M}\mathrm{A}}\right)\left(z_{\mathrm{F}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}\right)}{r}_{\mathrm{M}}>0.}}$(7)

Relatedness

We can calculate an individual’s consanguinity to any social partner (including themselves) as the probability that a gene drawn at random from the focal individual and a gene drawn at random from their social partner from the same locus (with replacement, in the event that the social partner is the focal individual themselves) are identical by descent (Bulmer, 1994). See Supporting Information for details of these calculations.

Religiosity optima

Setting the left-hand side of expressions (4) and (5) to zero and simultaneously solving yields the optimal level of religiosity for female and male juveniles as (8) $z_{\mathrm{F}\mathrm{J}}^{\ast }=z_{\mathrm{M}\mathrm{J}}^{\ast }={\displaystyle \frac{1+r_{\mathrm{s}\mathrm{i}\mathrm{b}}}{2}}$(8) if the ϵ_FJ and ϵ_MJ terms are neglected. Consideration of the ϵ_FJ and ϵ_MJ terms reveals that if both sexes adopt this optimal level of religiosity, then any perturbation from this value by either sex will result in selection acting to neutralize this perturbation, such that this optimum is stable. Accordingly, juveniles of both sexes are favored to exhibit the same level of religiosity, and hereafter we consider that this represents a single, non-sex-specific trait. Consequently, the optimal level of religiosity from the individual's perspective is z_{iJ |I}* = (1 + r_{sib |I})/2, or:$$ (9) $z_{\mathrm{i}\mathrm{J}|\mathrm{I}}^{\ast }={\displaystyle \frac{2d_{\mathrm{M}}\left(5-\alpha \right)\left(2-\text{β}\right)-d_{\mathrm{M}}^2\left(5-\alpha \right)\left(1-\text{β}\right)-d_{\mathrm{F}}^2\left(1-\alpha \right)\left(5+\text{β}\right)+8\left(\alpha +\text{β}\right)+2d_{\mathrm{F}}\left(10-d_{\mathrm{M}}\left(5-\alpha \right)+\text{β}-\alpha \left(6+\text{β}\right)\right)}{8\left(\alpha -d_{\mathrm{F}}^2\left(1-\alpha \right)+2d_{\mathrm{F}}\left(2-d_{\mathrm{M}}-\alpha \right)+2d_{\mathrm{M}}\left(2-\text{β}\right)-d_{\mathrm{M}}^2\left(1-\text{β}\right)+\text{β}\right)}.}$(9) From the maternal-origin gene's perspective the optimum is z_{iJ |M}* = (1 + r_{sib |M})/2, or: (10) $z_{\mathrm{i}\mathrm{J}|\mathrm{M}}^{\ast }={\displaystyle \frac{-3d_{\mathrm{F}}^2\left(1-\alpha \right)-3d_{\mathrm{M}}^2\left(1-\text{β}\right)+4\left(\alpha +\text{β}\right)+d_{\mathrm{M}}\left(12-\alpha -7\text{β}\right)+d_{\mathrm{F}}\left(12-7\alpha -\text{β}-d_{\mathrm{M}}\left(6-\alpha -\text{β}\right)\right)}{4\left(\alpha -d_{\mathrm{F}}^2\left(1-\alpha \right)+2d_{\mathrm{F}}\left(2-d_{\mathrm{M}}-\alpha \right)+2d_{\mathrm{M}}\left(2-\text{β}\right)-d_{\mathrm{M}}^2\left(1-\text{β}\right)+\text{β}\right)}.}$(10) From the paternal-origin gene’ s perspective the optimum is z_{iJ |P}* = (1 + r_{sib |P})/2, or: (11) $z_{\mathrm{i}\mathrm{J}|\mathrm{P}}^{\ast }={\displaystyle \frac{-d_{\mathrm{M}}^2\left(2-\alpha \right)\left(1-\text{β}\right)-d_{\mathrm{F}}^2\left(1-\alpha \right)\left(2+\text{β}\right)+4\left(\alpha +\text{β}\right)+d_{\mathrm{F}}\left(8-5\alpha -d_{\mathrm{M}}\left(4-\alpha +\text{β}\right)+3\text{β}-2\alpha \text{β}\right)+d_{\mathrm{M}}\left(8-3\alpha -3\text{β}+2\alpha \text{β}\right)}{4\left(\alpha -d_{\mathrm{F}}^2\left(1-\alpha \right)+2d_{\mathrm{F}}\left(2-d_{\mathrm{M}}-\alpha \right)+2d_{\mathrm{M}}\left(2-\text{β}\right)-d_{\mathrm{M}}^2\left(1-\text{β}\right)+\text{β}\right)}.}$(11) Setting the left-hand side of expressions (6) and (7) to zero and simultaneously solving yields the optimal level of religiosity for women as (12) $z_{\mathrm{F}\mathrm{A}}^{\ast }={\displaystyle \frac{\left(1-2z_{\mathrm{M}\mathrm{A}}+z_{\mathrm{M}\mathrm{A}}{r}_{\mathrm{F}}\right)+2r_{\mathrm{F}}+\sqrt{4\left(2+r_{\mathrm{F}}\right)z_{\mathrm{M}\mathrm{A}}+{\left(1-2z_{\mathrm{M}\mathrm{A}}+r_{\mathrm{F}}\left(2+z_{\mathrm{M}\mathrm{A}}\right)\right)}^2}}{2\left(2+r_{\mathrm{F}}\right)},}$(12) and for men as (13) $z_{\mathrm{M}\mathrm{A}}^{\ast }={\displaystyle \frac{\left(1-2z_{\mathrm{F}\mathrm{A}}+z_{\mathrm{F}\mathrm{A}}{r}_{\mathrm{M}}\right)2r_{\mathrm{M}}+\sqrt{4\left(2+r_{\mathrm{M}}\right)z_{\mathrm{F}\mathrm{A}}+{\left(1-2z_{\mathrm{F}\mathrm{A}}+r_{\mathrm{M}}\left(2+z_{\mathrm{F}\mathrm{A}}\right)\right)}^2}}{2\left(2+r_{\mathrm{M}}\right)}.}$(13)

Here again, optimal level of religiosity for adults of both sexes depends on relatedness to one’s parents’ social partners, expressed relative to relatedness to one’s parents, as well as the level of religiosity exhibited by adults of the opposite sex. As above, we can calculate the optima from the individual’s perspective (z_{iA |I}*), the maternal-origin gene's perspective (z_{iA |M}*) and the paternal-origin gene’s perspective (z_{iA |P}*) in adults by inserting the respective consanguinities. The levels of religiosity exhibited when both sexes are behaving optimally may be found by evaluating Equation (12) at z_MA = z_MA* and evaluating Equation (13) at z_FA = z_FA* and simultaneously solving for z_FA* and z_MA*; A numerical illustration of these optima is provided in Figures 1 and 2.