Genetic, maternal, and environmental influences on sociality in a pedigreed primate population
Study subjectsSubjects in our study are individually recognized wild capuchins found in and around the Lomas Barbudal Biological Reserve in Guanacaste, Costa Rica. This population has been under observation since 1990 (Perry 2012; Perry et al. 2012), including near continuous observation from January 2002 through March 2020.Data collectionWe use proximity data on subjects collected during group scan sampling between January 2001 and March 2020 (Altmann 1974). Included in scans are the identity of the subject, and the identity of other individuals within approximately 4 meters of them. Scans have been collected on all individuals in study groups since 2002, and on all adults and subadults since 2001. Scans are taken opportunistically, without regard to time of day. At least 10 min separate consecutive scans of the same individual to reduce the non-independence of scans taken close in time.Data in this manuscript were collected by 124 observers, with an average of 7.1 data collectors per month. Observers typically work in teams of two to three and rotate across different groups to reduce potential observer bias. Observers also rotate across observer teams to avoid observer drift in coding, since observer teams could potentially start to code behaviors differently from each other in the absence of overlap in observer composition.Initial pedigree constructionOf the 376 individuals in our behavioral dataset, 280 (74.5%) were first seen within three months of their births, and we could confidently assign maternity to them based on demographic (pregnancies) and behavioral data (primary nursing) even prior to genotyping. Of the remaining individuals, 41 (10.9%) were males of unknown origin that immigrated into our study population, while the rest were natal to our study groups but were first seen as older infants ( >3 months), juveniles, or (sub)adults (14.6%) and required genotyping to assign/confirm maternity. Paternity was assigned based on genetic information when possible (but see Non-genotyped individuals).In total, 287 subjects (76.3%) had two assigned parents, 37 had one assigned parent (9.8%), and 52 (13.8%) had no assigned parent based on demographic, behavioral, and/or genetic parentage information. Most individuals with no assigned parents were immigrant males (78.9%).GenotypingInformation on genetic parentage assignment (at up to 18 microsatellite loci) in our study population is available from previously published work (1996–2005 (Muniz et al. 2006), 2005–2012 (Godoy et al. 2016b)). Partial genotypes (up to 14 loci) have been generated for individuals in this study which more recently entered the study population through birth or immigration (n = 91, 2012–2020) (See SI File 1). Briefly, DNA was extracted primarily from non-invasively collected fecal samples, and occasionally from tissue samples obtained from deceased individuals, then amplified at up to 18 autosomal tetranucleotide microsatellite loci (Muniz and Vigilant 2008) using either a 1-step or 2-step PCR protocol (Arandjelovic et al. 2009). There were no significant deviations from Hardy-Weinberg equilibrium, and no evidence of linkage disequilibrium between loci was found (Muniz 2008).DNA samples were run at a minimum in triplicate, but additional PCRs were performed on low quality samples (e.g., with low quantities of DNA). Genotypes at each of the loci were assigned to be heterozygous when each allele was seen at least twice in independent PCRs, and assigned as homozygous when the allele was seen in at least three independent PCRs in absence of a second allele.Amplicons were analyzed using an ABI PRISM3100 automated sequencer and GeneMapper Software (Applied Biosystems, Foster City, CA, USA). Likelihood-based parentage assignments were performed using CERVUS 2.0 or 3.0 (Marshall et al. 1998; Kalinowski et al. 2007). The average exclusionary power of the 18 microsatellites was 0.9888 for the first parent and 0.9998 for second parent (Muniz et al. 2006).Individuals with unknown parents (e.g., immigrant males, founders) were genotyped twice (i.e., using two independent DNA samples) following the procedures described above to guard against sample mix up. Known mother-offspring pairs were confirmed by ascertaining the absence of Mendelian mismatches across all loci for the pair, though one mismatch was allowed to account for null alleles, mutations, and genotyping errors. We detected one null allele in the population in 19 individuals and traced it back to a male who was either the father or grandfather of those individuals (Muniz et al. 2006; Godoy et al. 2016b).Candidate males for paternity assignment were chosen based on group membership around the time of an infant’s conception (typically 1–10 males). In cases when conceptions occurred prior to the habituation of a study group, we used the identities of all adult males present when the group was first observed. Candidate mothers were similarly chosen for individuals that were first seen as older infants, juveniles, or (sub)adults. For individuals born post-group habituation, CERVUS has always assigned paternity from the pool of potential candidate fathers. Parent-offspring pairs and trios were allowed one mismatch (excluding those at the locus with the known null allele).Pedigree updatingNon-genotyped individualsDuring stable tenures, alpha males in our population sire approximately 73% of infants born in their groups, including 90% of offspring born to unrelated females (Godoy et al. 2016a). There is strong evidence of inbreeding avoidance between alpha males and their female descendants, with relatedness to females as the primary factor constraining alpha male monopolization of paternity within groups (Muniz et al. 2006, 2010; Godoy et al. 2016a, 2016b; Wikberg et al. 2017, 2018). We used this information to update our pedigree, filling missing father information with the identity of the alpha male around the time of a non-genotyped individual’s conception, but only if their mother was not the daughter or granddaughter of the alpha male (i.e., with inbreeding avoidance). This approach allowed us to assign presumed paternity to 21 non-genotyped individuals (5.6% of subjects) who were natal to our study groups.Individuals with missing or incomplete parentageOut of the original four study groups (from which fissions led to eight additional study groups), we lacked parentage information (i.e., neither parent was sampled) for 12 individuals first seen at the time of habituation. We had incomplete parentage on an additional 11 adults (i.e., only one parent was sampled). We used the software program COLONY version 2.0.6.7 to look for evidence of whether these individuals were related to each other at the level of full sibling (Jones and Wang 2010). We also looked for potential full sibling pairs among the non-natal immigrant males in the population, since co-migrant males are typically kin (Perry 2012; Wikberg et al. 2014, 2018). We assigned five full sibling pairs among co-migrant males, and four full sibling pairs among natal founders. For any remaining co-migrant males and natal founder pairs that were not assigned as full siblings, we assumed these to be either paternal (migrants) or maternal (natal) half siblings, as is typical in this study population (Perry 2012). These assignments are likely to have some error. However, based on what we know about kinship in capuchins, it would introduce more error to assume that these pairs are unrelated.We pruned our modified population pedigree using the R package pedantics version 1.01 (Morrissey and Wilson 2010), to include only individuals that were linked to the subjects in our behavioral dataset. The reduction in missing data can improve convergence and mixing of models (Hadfield 2010). The pruned pedigree contained 419 individuals, with 353 maternities, 354 paternities, 209 full sibships, 413 maternal half sibships, and 1496 paternal half sibships. Maximum pedigree depth was six generations (mean = 3.03).Sociality measures (response variables)We generated two related proximity-based measures of sociality—(1) whether an individual was seen in proximity of another monkey (within ~4 meters) during a scan (i.e., they were not alone), and (2) the number of partners an individual has nearby (within ~4 meters) during a scan. The former is measure of the propensity of an individual to be social versus alone, while the latter is more indicative of the gregariousness of an individual. These two phenotypes are not independent, as they are generated from the same data (Fig. 1a).Fig. 1: Distribution of sociality, sampling, group size, and alpha tenure length.The scatterplot in a shows the proportion of scans per individual per month where the subjects were recorded in proximity of others on the x-axis, and the average number of social partners per scan per month for subjects on the y-axis. The sizes of the circles in a are proportional to sample size (range: 5–317 scans per data point). The figure in b shows the number of calendar years of data sampling per subject (range: 1–20), c variation in group size, and d the number of calendar years represented by different alpha tenures in the dataset. Note that d does not represent the full diversity of alpha tenure lengths in the population, only within the dataset: some tenure lengths are left-truncated as data from 1990–2000 are not included in this dataset. Figure produced in R using ggplot2 version 3.3.5 (Wickham 2016) and cowplot version 1.1.1 (Wilke 2020). The capuchin image was generated in R using sketcher version 0.1.3 (Tsuda 2020) based on an image taken by Nicholas Schleissmann.Full size imageWe compiled the scans of individuals by month (mean: 31.9, range: 5–317 scans per month) so that we had counts of (1) the total number of scans where an individual was social and (2) the total number of partners an individual had. With these counts we could look at the (1) proportion of time spent social (versus alone) and (2) the average number of partners an individual had, while still preserving information about sampling density (number of scans).To be included in any month, subjects needed to have at least five scan samples in that period. As we are interested in the repeatability of our measures of social behavior, subjects had to have at least six months of data to be included.We excluded dependent infants (less than one year of age) as potential social partners of their mothers. We also excluded these dependent infants as subjects, since an infant is expected to be in close proximity of its mother, particularly during the first half of their first year of life (Godoy 2010; Perry 2012). Including data from infants would likely introduce upward bias to heritability estimates, because mothers and their dependent offspring (whom share high relatedness) would often be in close proximity of each other, and their measures of proximity to others would thus also be highly correlated.On average, subjects spent just over half of their sampled time within approximately four meters of another monkey (mean: 0.539, standard deviation: 0.193) and had approximately one social partner per scan (mean: 1.057, standard deviation: 0.619) (Fig. 1a). Our dataset consisted of 22,138 monthly sociality scores on 376 subjects generated from 641 140 scans (mean: 56.5 months per subject, range: 6–184 months per subject). Almost all subjects (99.7%, i.e., all but one) were represented by data across more than one calendar year (25, 50, 75% quantiles: 4, 7, 10 different years of data collection) (Fig. 1b).Fixed effectsWe included age (as a cubic function) and sex in our models, as well as their interaction to account for differences in how male and female capuchins sexually mature and age. Age in our dataset was right-skewed with higher representation at younger ages (mean: 9.3 years, standard deviation: 6.9) (Fig. 2). To put the ages in developmental context, mean age at first live birth is around 6.3 years for females in this population, though females can begin reproducing in their 5th year (Perry et al. 2012). Males as young as six years old have been known to sire offspring (Godoy et al. 2016b), but males tend to not reach full adult size until their 10th year (Jack et al. 2014).Fig. 2: Sociality as a function of age and sex.Circles represent individual monthly data. The sizes of the circles are proportional to sample size (range: 5–317 shows per data point). Circles in a represent the proportion of time individuals were seen in proximity of others (not alone) per month, while in b represent the average number of partners for individuals per month. Solid lines represent estimated sociality scores based on age and sex, with all other fixed effects set to the mean. The two x-axes represent age as z-scores and in years. Figure produced in R using ggplot2 version 3.3.5 (Wickham 2016).Full size imageSeasonal environmental changes, such as in food abundance, or temperature and rain, can lead to changes in how individuals cluster near others, for example, because of how food resources become distributed in the environment. For example, in black-crested gibbons (Nomascus concolor), group averages of dyadic proximity have been documented to decrease from the dry season to wet season, with increased average group proximity during cold months and lowered proximity during warm months (Guan et al. 2013). We account for seasonal variation by modelling monthly changes as a sine wave, through inclusion of the sine and cosine functions of a transformed month variable (See SI File 1 for further details).Central America is a region of ENSO-related precipitation, where the El Niño-Southern Oscillation (ENSO) has an impact on large scale patterns of temperature and precipitation (Ropelewski and Halpert 1987). Bimonthly rainfall anomalies are linked with both the warm El Niño and cool La Niña phases in a neighboring tropical dry forest in Costa Rica, where long-term monitoring of wild white-faced capuchins has shown declines in reproductive output associated with El Niño-like conditions (Campos et al. 2015). To account for the large-scale influence of ENSO on group dynamics, we included a factor variable for three different ENSO phases (Average/Neutral, Cool/La Niña, and Warm/El Niño). We used the bi-monthly Multivariate El Niño/Southern Oscillation (ENSO) index (MEI.v2) obtained from the Physical Sciences Laboratory of the National Oceanic and Atmospheric Administration (https://psl.noaa.gov/enso/mei/, retrieved: 2021-11-06) to determine the different phases. MEI.v2 is a composite index of five different variables (sea level pressure, sea surface temperature, surface zonal winds, surface meridional winds, and Outgoing Longwave Radiation) used to create a time series of ENSO conditions from 1979 to present (Zhang et al. 2019). Warm phases correspond to MEI.v2 values of 0.5 or higher, while cool phases correspond to values of −0.5 or lower.Demographic differences between groups and within groups across time can also lead to variation in behavior. For example, group size has been found to correlate with the amount of time that individuals spend grooming in various primate species (Dunbar 1991; Lehmann et al. 2007). Group size is also associated with higher sociality measures such as both the number of strong and weak ties that individuals form in diverse clades of primates (Schülke et al. 2022). We attempt to account for variation that arises from such demographic differences by including group size (mean: 24.7, standard deviation: 7.9) (Fig. 1c) as a fixed effect.In our models, group size and cubic age were centered and scaled to a mean of zero and a standard deviation of one.Random effectsAll models include the identity of the subject (VID, n = 376) as a random factor, as well as subject identity nested within year (VID:Year, n = 3150), the identity of each subject’s mother (VM, n = 142), maternal identity nested within group of residence within year of data collection (VM:GroupAlpha:Year, n = 2085), and a special variable known as the animal term to account for additive genetic variance (VA). These components contribute to long- and/or short-term repeatability of individuals. All models also include year of data collection (VYear, n = 20), month nested within year (VMonth:Year, n = 224), and the identity of each subject’s group of residence both across years (VGroupAlpha, n = 56) and within years (VGroupAlpha:Year, n = 200).VID in the models (since the models also additionally estimate VM and VA) can be thought of as estimating the “permanent environment variance” (i.e., VPE) of an individual, which is the “individual-specific variation in environmental conditions that permanently affect the phenotype (e.g. early-life conditions)” (Dingemanse et al. 2010). VID:Year captures the variance explained by the repeated sampling of the same individuals within a particular year. We use it to estimate the proportion of the phenotypic variance due to similarity in the trait within individuals from data taken closer in time (within the same year). During such a relatively short period, individuals are more likely to be stable in important social traits such as kin availability, dominance rank for adults, and maternal dominance rank for infants and young juveniles.VM estimates the variance explained by maternal effects (m2), specifically similarity between maternal siblings. Maternal identities were not available for all subjects, namely 11 immigrant males of unknown origin who were not assigned by COLONY as having a full sibling. We created unique dummy codes for their maternal identities, so that no two of these individuals shared the same mother. We additionally nested maternal identities (VM:GroupAlpha:Year) to account for similarity between maternal siblings residing in the same group in the same year. Such a nested structure might capture potential upward biases on heritability due to maternal kin biases in spatial association among siblings residing in the same group.We estimate h2 in our models by fitting a random effects term (VA), referred to as the animal term, which in the R package MCMCglmm links to the identities of individuals in our population pedigree (Hadfield 2010; see below for details on the implementation of the models in MCMCglmm). Inclusion of the animal term provides our models with an additive genetic variance component based on the estimated coefficients of relatedness between individuals in our pedigree. In short, if animals that share more alleles are also more like each other in their behavior, then variation in the behavior may well be due to genetic variation in the population (under the assumption that phenotypic similarity is not due to a shared environment, or is adequately controlled for by fixed and random effects in the model).VYear and VMonth:Year were included in order to account for temporal variation in sociality scores not captured by the fixed effects of seasonality or ENSO phase. These could arise from, for example, observer drift in coding (i.e., measurement error) or prevailing environmental conditions (e.g., drought) that could lead to changes to how individuals cluster near others. There were 218 unique observer combinations across the 224 months represented in the dataset, so VMonth:Year should also capture variance due to any differences between observer teams, though we cannot separate out the unique influence of observers.VGroupAlpha represents variance arising from the different alpha tenures within groups in our study population. VGroupAlpha captures both variance due to group of residence effects and the additional influence of alpha tenures within those groups. In capuchins, alpha males are ‘keystone’ individuals, whose influence is disproportionate relative to that of others in the population, and thus play important roles in establishing group dynamics (Jack and Fedigan 2018). Including group of residence, as defined by alpha tenure, is also important because it helps to account for the higher relatedness within groups within alpha tenures which results from high male reproductive skew toward alpha males. At Lomas Barbudal, males can remain in their alpha position for upwards of 18 years. Alpha tenures in this dataset spanned one to 14 years (Fig. 1d), so we additionally nested the identity of alpha males per group within years (VGroupAlpha:Year) so as to separate the within-year and across-year influences of group of residence.Statistical methodsWe ran analyses in R 4.1.2 (R Core Team 2021), using a Bayesian method with the R package MCMCglmm version 2.32 (Hadfield 2010). Data and code used to run all models is provided in the Supplementary Information.For our binary response variable (social versus alone), which was pooled into monthly units, we fit models with a binomial distribution and logit link function (family = “multinomial2”), with the number of scans each individual was documented social (‘successes’) versus the number of times alone (‘failures’).For our other response variable (number of partners), which was also pooled into monthly units, we fit models with a Poisson distribution (family = “poisson”), with the total number (sum) of partners per month. We included the natural log of the number of scans per month as a fixed effect to account for sampling effort. We set a strong prior for the log of sampling effort so that the rate at which events occurred was 1 (i.e., we could look at average number of partners per scan).We used a parameter-expanded prior (V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000) and two inverse Wishart priors (V = 1, nu = 0.002; V = 1, nu = 0.02) for the G structures in our models (i.e., random effects variance components). The prior on the residual variance component was set to one for both the binomial and Poisson models. Estimates for variance components were robust against the choice of prior (SI Fig. 3). We therefore only report findings from models run with parameter-expanded priors in the main text.Pilot runs (thin = 10, burnin = 3000, nitt = 13,000) indicated that autocorrelation values would remain high for some variance components in models run with parameter-expanded priors, even with large thinning intervals. We therefore increased the number of iterations to guarantee effective sample sizes of at least 1000, but ideally closer to 4000. All models were run with a long burn-in period of at least 10,000 iterations.We ran multiple chains (n = 4) of each model and assessed convergence of the chains visually (SI Files 2a-b), as well as through the Gelman-Rubin criterion implemented via the ‘gelman.diag’ function from the coda package in R (version 0.19-4) (Plummer et al. 2006). Scale reduction factors were below 1.02, signifying good convergence. We used Heidelberger and Welch’s convergence diagnostic test for stationarity to check convergence of each chain using the ‘heidel.diag’ function from the coda package. Results are presented from the first chain of each model.Reduced modelsInclusion of fixed effects can potentially have an impact on the estimates of variance components in models because total phenotypic variance (VP) is estimated (and partitioned among the different random effects) after conditioning on the fixed effects. Heritability estimates, for example, can be higher because the variance explained by the fixed effects structure (VFE) is not included in VP, thus making the relative contribution of VA to VP larger compared to the same model without fixed effects (Wilson 2008). Conversely, not adequately controlling for relevant fixed effects that contribute to phenotypic variance among and within individuals may potentially lead to an underestimation of VA and associated heritability (h2).We ran multiple reduced versions of our models to look at the impact of fixed effects on our variance components. We began with an intercept-only version (i.e., no fixed effects), then built-up complexity by adding in versions with the properties of the individuals first (age, sex), then properties of the group (group size), and subsequently environmental properties (seasonality, ENSO phases). Outputs for these reduced models are provided in the Supplementary Information (SI Table 2, SI Table 3).We provide the deviance information criterion (DIC) values for models (automatically generated by the MCMCglmm package). DIC is a generalization for multi-level models of the Akaike Information Criterion (AIC); and as in AIC, lower DIC values indicate better fit.Transformations from unobserved latent scale to observed data scaleOutputs from our MCMCglmm models were on the unobserved latent scale. We used the R package QGglmm (version 0.7.4) to additionally compute parameters of interest on the observed data scale (de Villemereuil et al. 2016; de Villemereuil 2018). We used the functions ‘QCicc’ to compute Intra-Class Correlation (ICC) coefficients and ‘QGparams’ to compute additive genetic variance and thus narrow-sense heritability (h2) on the observed data scale. We implemented the ‘QGparams’ and ‘QGicc’ functions with parameters model = ‘binomN.logit’ and n.obs = 32 (the average number of scans per subject per month in our dataset) for the binomial model and model = ‘Poisson.log’ for the Poisson model. The choice of value for n.obs is somewhat arbitrary, and we show the consequences for changes in values of this parameter (i.e., higher estimates with increasing values of n.obs) in SI Fig. 4.Closed form solutions in QGglmm are not available for integrating over posterior distributions generated from binomial models with logit link functions (de Villemereuil 2016). Consequently, using the ‘QGicc’ function is particularly slow. We therefore estimate ICCs from our binomial models using a random subset of the posterior (n = 1000 iterations).The code used for transforming the MCMCglmm outputs from the latent scale to the original data scale are available online (see DATA AVAILABILITY).Repeatability and the proportion of variance explained by variance componentsTotal phenotypic variance (VP) was the sum of estimates from all variance components and residual variance in a model (VP = VID + VID:Year + VM + VM:GroupAlpha:Year + VA + VGroupAlpha + VGroupAlpha:Year + VMonth:Year + VYear + Vresidual). The proportion of variance explained by each variance component was calculated by including its estimate in the numerator while including total phenotypic variance in the denominator. So, for example the proportion of variance explained by year of data collection was calculated as (left( {frac{{V_{Year}}}{{V_P}}} right)).Long-term repeatability was calculated with the sum of VID, VM, and VA in the numerator. Short-term repeatability was calculated similarly but with inclusion of within-series variances (VID + VM + VA + VID:Year + VM:GroupAlpha:Year) in the numerator to capture additional consistency in among-individual differences resulting from greater environmental similarity within a time series (i.e., year).We report posterior modes and 95% Highest Posterior Density intervals (i.e., 95HPDI in square brackets). Unless mentioned otherwise, we present results on the unobserved latent scale, and without the variance from the fixed effects (VFE) incorporated into VP. For completeness, estimates with VFE included in VP and transformations to the observed data scale are also provided in SI Table 3. More
