The genomics of ecological flexibility, large brains, and long lives in capuchin monkeys revealed with fecalFACS.

Ecological flexibility, extended lifespans, and large brains have long intrigued evolutionary biologists, and comparative genomics offers an efficient and effective tool for generating new insights into the evolution of such traits. Studies of capuchin monkeys are particularly well situated to shed light on the selective pressures and genetic underpinnings of local adaptation to diverse habitats, longevity, and brain development. Distributed widely across Central and South America, they are inventive and extractive foragers, known for their sensorimotor intelligence. Capuchins have among the largest relative brain size of any monkey and a lifespan that exceeds 50 y, despite their small (3 to 5 kg) body size. We assemble and annotate a de novo reference genome for Cebus imitator Through high-depth sequencing of DNA derived from blood, various tissues, and feces via fluorescence-activated cell sorting (fecalFACS) to isolate monkey epithelial cells, we compared genomes of capuchin populations from tropical dry forests and lowland rainforests and identified population divergence in genes involved in water balance, kidney function, and metabolism. Through a comparative genomics approach spanning a wide diversity of mammals, we identified genes under positive selection associated with longevity and brain development. Additionally, we provide a technological advancement in the use of noninvasive genomics for studies of free-ranging mammals. Our intra- and interspecific comparative study of capuchin genomics provides insights into processes underlying local adaptation to diverse and physiologically challenging environments, as well as the molecular basis of brain evolution and longevity.

Large brains, long lifespans, extended juvenescence, tool use, and problem solving are hallmark characteristics of great apes, and are of enduring interest in studies of human evolution (1–4). Similar suites of traits have arisen in other lineages, including some cetaceans, corvids and, independently, in another radiation of primates, the capuchin monkeys. Like great apes, they have diverse diets, consume and seek out high-energy resources, engage in complex extractive foraging techniques (5, 6) to consume difficult-to-access invertebrates and nuts (6), and have an extended lifespan, presently recorded up to 54 y in captivity (7, 8). While they do not show evidence of some traits linked with large brain size in humans (e.g., human-like social networks and cultural and technological transmission from older to younger groupmates), their propensity for tool use and their ecological flexibility may have contributed to their convergence with the great apes (9), offering opportunities for understanding the evolution of key traits via the comparative method (10–12). Similar approaches have revealed positive selection on genes related to brain size and long lives in great apes and other mammals (13, 14), but our understanding of the genetic underpinnings of these traits remains far from complete.

Capuchins also offer excellent opportunities to study local adaptation to challenging seasonal biomes. They occupy diverse habitats, including rainforests and, in the northern extent of their range, tropical dry forests. Particular challenges of the tropical dry forest are staying hydrated during the seasonally prominent droughts, high temperatures in the absence of foliage, and coping metabolically with periods of fruit dearth (Fig. 1). The long-term study of white-faced capuchins (Cebus imitator) occupying these seasonal forests has demonstrated that high infant mortality rates accompany periods of intense drought, illustrating the strength of this selective pressure (15). Furthermore, the seasonally low abundance of fruit is associated with muscular wasting and low circulating levels of urinary creatinine among these capuchins (16). Additionally, the sensory challenges of food search in dry versus humid biomes are also distinct. Odor detection and propagation is affected by temperature and humidity (17), and color vision is hypothesized to be adaptive in the search for ripe fruits and young reddish leaves against a background of thick, mature foliage (18), which is absent for long stretches in dry deciduous forests. The behavioral plasticity of capuchins is widely acknowledged as a source of their ability to adapt to these dramatically different habitats (19–21). However, physiological processes, including water balance and metabolic adaptations to low caloric intake, and sensory adaptations to food search, are also anticipated to be targets of natural selection, as seen in other mammals (22–24). Understanding population-level differences between primates inhabiting different biomes, contextualized by their demographic history, genomic diversity, and historical patterns of migration, will generate new insights.

Fig. 1.

SSR during wet (Left) and dry (Center) seasons. (Right) Map of sampling locations in Costa Rica. The two northern sites, SSR and Cañas, have tropical dry-forest biomes, whereas the two southern sites, Quepos and Manuel Antonio, are tropical wet forests. Photos courtesy of A.D.M. Drawing of white-faced capuchin monkey by Alejandra Tejada-Martinez; map courtesy of Eric Gaba–Wikimedia Commons user: Sting.

Unfortunately, high-quality biological specimens from wild capuchins are not readily available. As is the case with most of the world’s primates, many of which are rare or threatened (25), this has limited the scope of questions about their biology that can be answered. Although recent advances in noninvasive genomics have allowed for the sequencing of partial genomes by enriching the proportion of endogenous DNA in feces (26–29), it has not yet been feasible to sequence whole genomes from noninvasive samples at high coverage; this has limited the extent to which noninvasive samples can be used to generate genomic resources for nonmodel organisms, such as capuchins.

Toward identifying the genetic underpinnings of local adaptation to seasonally harsh environments, large brains, and long lifespans, we assembled and annotated a reference genome of C. imitator (). Additionally, we sequenced the genomes of individuals inhabiting two distinct environments in Costa Rica: Lowland evergreen rainforest (southern population) and lowland tropical dry forest (northern population). We conducted high-coverage resequencing (10× to 47×) for 10 of these individuals, and sequenced an additional 13 at low-coverage (0.1× to 4.4×). Importantly, to facilitate the population-wide analyses without the need for potentially harmful invasive sampling of wild primates, we developed a method for minimally biased, whole-genome sequencing of fecal DNA using fluorescence-activated cell sorting (fecalFACS) that we used to generate both high- and low-coverage genomes (Fig. 2). With these genomes, we assess the genetic underpinnings of capuchin-specific biology and adaptation in a comparative framework. First, we scanned the high-coverage genomes (six from the northern dry forest and four from the southern rainforest) for regions exhibiting population specific divergence to assess the extent of local adaptation to dry forest and rainforest environments. We examine how genes related to water balance, metabolism, muscular wasting, and chemosensation have diverged between populations. Second, we conduct an analysis of positive selection on the white-faced capuchin genome through codon-based models of evolution and enrichment tests focusing on genes that may underlie brain development and lifespan. Third, we identify the population structure, genomic diversity, and demographic history of the species using a mixture of traditional and noninvasive fecalFACS genomes (n = 23).

Fig. 2.

Mapping percentages of sequencing reads from RNAlater preserved fecal DNA libraries prepared with FACS for (A) all samples (box-plot elements: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers), and (B) individual libraries. (C) Increase in mapping rate for RNAlater preserved samples. (D) Relationship between mapped read duplication and number of cells with LOESS smoothing. The duplication rate decreases sharply once a threshold of about 1,000 cells is reached.

Results

Local Adaptation to Seasonal Food and Water Scarcity.

We predicted that genes related to water balance, metabolism, muscular wasting, and chemosensation would differ between dry forest and rainforest populations of white-faced capuchins, reflecting local adaptation to different habitats. To test this, we searched for associations between genes in windows with high FST and divergent nonsynonymous single nucleotide variants (SNV) with high or moderate predicted effects. Of the 299 genes identified in high FST windows, 39 had highly differentiated (FST ≥ 0.75) nonsynonymous SNVs ( and Dataset S1). We identified 26 genes with nonsynonymous SNVs of high or moderate effect that are fixed between populations in our dataset; 13 of these overlapped with the set of 39. Our enrichment analysis identified a single significant gene ontology (GO) biological process: Regulation of protein localization to cilium (GO: 1903564). In accordance with our hypothesis, disruptions of cilia proteins are predominantly associated with disorders of the kidney and retina (ciliopathies) (30). Furthermore, enrichment of genes associated with disease states was linked to kidney function, metabolism, and muscular wasting. These genes are good candidates for adaptive resilience to seasonal water and food shortages and warrant further investigation. We highlight several genes of particular promise (Fig. 3 and ).

Fig. 3.

Highly differentiated genes between wet- and dry-forest populations involved in diabetes, kidney function, and creatinine levels. (A) Hudson’s FST within windows of 20 kb with a 4-kb slide. Gene regions are in red, flanked by 500 kb (or length to beginning or end of scaffold) of sequence. The x axis values correspond to position along the scaffold. The dotted line indicates the average FST value across all windows (FST = 0.178), and the dashed line represents the top 0.5% of values (FST = 0.797). Vertical black lines indicate a nonsynonymous SNP with an Fst ≥ 0.750, excluding BCAS3 (see ). (B) Heatmaps indicating the pattern of SNP variation within and surrounding highly divergent genes. SNVs within the genes are located under the red band and those within 200 kb of flanking region under the gray bands.

Evidence of adaptation to food and water scarcity.

Population differences in multiple candidate genes indicate that dry-forest capuchins could be adapted to seasonal drought-like conditions and food scarcity. SERPINC1 encodes antithrombin III, which is involved in anticoagulant and antiinflammatory responses associated with numerous kidney-related disorders including salt-sensitive hypertension, proteinuria, and nephrotic syndrome (31–34). Additionally, sequence variants of AXDND1 (as identified in the GeneCards database) are associated with nephrotic syndrome. BCAS3 is expressed in multiple distal nephron cells types (35), and is associated with four pleiotropic kidney functions (concentrations of serum creatinine, blood urea nitrogen, uric acid, and the estimated glomerular filtration rate based on serum creatinine level) (36). Although we did not identify any population-specific nonsynonymous SNVs in BCAS3, the genomic windows encompassing the gene rank among the highest regions of FST in our dataset, and intronic variation in BCAS3 putatively impacts estimated glomerular filtration rate in humans (35).

We also examined candidate genes associated with metabolism and with catalysis of muscle tissues, which occurs when energy intake is insufficient to fuel metabolic and energetic needs over long periods. We found population differentiation in genes associated with congenital muscular dystrophies—including mild forms that present with muscular wasting—and abnormal circulating creatine kinase concentration (HP: 0040081): ITGA7, ISPD (CRPPA), and SYNE2 (37–40). The gene LAMA2 (41), which falls in a high FST window, has also been associated with muscular pathologies, some of which can be reduced by transgenic overexpression of ITGA7. Finally, we examined genes with roles in carbohydrate and lipid metabolism. We observed differentiation in GLIS3, which is associated with type 1, type 2, and neonatal diabetes (42). IRF9, which also has a fixed SNV of moderate effect between populations, has been shown to influence glucose metabolism, insulin sensitivity, and obesity in knockout mice (43). IRF9 has also been associated with immunity and susceptibility to viral infections (OMIM: 618648).

Given the appearance of both diabetes and kidney disorders in our gene sets, we conducted an a posteriori search of our high FST gene set for overrepresentation of genes associated with diabetic nephropathy (EFO_0000401) in GWAS Catalog (https://www.ebi.ac.uk/gwas) (). Seven genes were present in both our gene set of 299 high FST genes and the diabetic nephropathy set of 117. Given the 16,553 annotated genes with Hugo Gene Nomenclature Committee (HGNC) Gene IDs, 7 overlapping genes would occur with P = 0.00046 when permuted 100,000 times.

Evidence of adaptation in sensory systems.

Intact gene sequences were found for two opsin genes, OPN1SW and OPN1LW, which underlie color vision and acuity. In both populations, we observed similar levels of polymorphism in the long-wavelength sensitive opsin gene (OPN1LW), with the same allelic variants at each of the three long wave cone opsin tuning sites (180, A/S; 277, F/Y; and 285, T/A) (). We observed some evidence for population-specific variation associated with the photoreceptive layers of the retina; a fixed nonsynonymous SNV in CCDC66, which falls in a high FST region, is highly expressed in photoreceptive layers of the retina (44). Turning to olfaction, we identified 614 olfactory receptor (OR) genes and pseudogenes in the capuchin reference genome: 408 intact, 45 truncated, and 161 pseudogenized (). To test for population differences in the OR gene repertoire, we assembled each olfactory gene/pseudogene independently in each individual. The proportion of total functional ORs was stable across individuals and populations, with seemingly trivial fluctuations (North x̅ = 411, s = 1.6; South x̅ = 408.5, s = 1.3), possibly driven by a small difference in OR family 5/8/9 (). We also identified 7 intact vomeronasal (VR) and 28 taste receptor (TASR) and taste receptor-like genes (), two of which (TAS1R and TAS2R4) have nonsynonymous SNVs with fixed variants in the north ().

Evolution of Brain Development and Longevity in Capuchins.

To identify genes that underlie brain development and lifespan, which are of particular interest given the derived features of capuchin biology, we ran a codon-based positive selection analysis with PAML and subsequently tested for functional enrichment with ToppFunn (45) (Fig. 4). Of the genes identified as being under positive selection in the branch-site model (746 genes) and the branch model (612 genes), we identified significant enrichment in genes associated with brain development (52 genes) (Dataset S1), neurogenesis (77 genes) (Dataset S1), and microcephaly (31 genes) (Dataset S1). For example, WDR62, BPTF, BBS7, and NUP133 mutations are directly associated with brain size and related malformations, including microcephaly (46–49). Mechanistic target of rapamycin (MTOR) signaling malfunction is also implicated in developmental brain malformations (50), and NUP133 is involved in nuclear migration during mammalian brain development (51). Several genes are also linked with cognitive ability [e.g., PHF8 (52)].

Fig. 4.

Genes under positive selection in white-faced capuchin monkeys show enrichment associated with longevity (maximum recorded lifespan in captivity) and brain size and development. Several genes showing evidence of positive selection in the Cebus lineage are listed for each trait. The displayed species are the primates used for the PAML run, or a congeneric species in cases of missing trait data (e.g., C. capucinus in place of C. imitator). Relative brain size, calculated as EQ = brain mass/(0.085 × (body mass0.775)) (1), is displayed to account for the large range in body mass. Trait data are from refs. 1, 8, 109, and 110. Cebus capucinus silhouette credit: Phylopic/Sarah Werning, licensed under CC BY 3.0. Nomascus sp. silhouette credit: Phylopic/Kai R. Caspar, licensed under CC BY 3.0. Microcebus murinus silhouette credit: Phylopic/Maky, Gabriella Skollar, and Rebecca Lewis.

We found 48 genes putatively linked to longevity, as identified in the GenAge and CellAge databases (14, 53), to be under positive selection in capuchins, including PARP1, MTOR, SREBF1, INSR, HTT, RB1, and MDM2 (). Of note, poly (ADP ribose) polymerase 1 (PARP1) putatively serves as a determinant of mammalian aging due to its activity in the recovery of cells from DNA damage, and MTOR acts as a regulator of cell growth and proliferation while also being generally involved in metabolism and other processes. Additional key genes in aging and metabolism include sterol regulatory element binding transcription factor 1 (SREBF1), which acts as a regulator of the metabolic benefits of caloric restriction (54, 55), and the insulin receptor (INSR), a major player in longevity regulation (56). As for specific age-related diseases, huntingtin (HTT) is under selection in mammals; HTT is not only involved in Huntington’s disease but has also been associated with longevity in mice (57). Finally, various cell cycle regulators (e.g., RB1, MDM2) are also under positive selection in capuchins, and indeed, cell cycle is an enriched term among positively selected genes (), though these could be related to other life-history traits, like developmental schedules that correlate with longevity.

Population Genomics of Costa Rican Capuchins with FecalFACS.

The pattern of clustering in our maximum-likelihood SNV tree recapitulates the expected patterns of geographic distance and ecological separation in our samples (Fig. 5). Similarly, in the projected principal component analysis (PCA), all individuals from the seasonal dry forests in the northwest are sharply discriminated from individuals inhabiting the southern rainforests along PC1. Levels of heterozygosity calculated in overlapping 1-Mb genomic windows (with a step size of 100 kb) were significantly higher in the southern population (W = 1,535,400,000, P < 2.2e-16) (Fig. 6 and ). Furthermore, the median pairwise heterozygosity for each southern individual (range: 0.00065 to 0.00071) was higher than any northern monkey (0.00047 to 0.00057) (W = 0, P = 0.009524) (). In the northern population, we also identified long runs of homozygosity significantly more often (W = 24, P = 0.009524), and more of the longest runs (≥5 Mb) (W = 1315.5, P = 0.03053) (Fig. 6 and ). Pairwise sequential Markovian coalescent (PSMC) analysis of demographic history (Fig. 6) reveals that white-faced capuchins had a peak effective population size of ∼60,000 effective individuals ∼1 Ma, which declined to fewer than 20,000 during the middle to late Pleistocene transition. After recovering during the middle Pleistocene, they declined precipitously through the late Pleistocene to fewer than 5,000 effective individuals.

Fig. 5.

Population subdivision in C. imitator. (Left) PCs of 13 fecal and 10 blood/tissue libraries from white-faced capuchins. Individuals from northern and southern sites separate on PC1. Low- and high-coverage C. imitator samples from Santa Rosa plot in the same cluster. (Right) Maximum-likelihood tree of 9 fecal and 10 blood/tissue libraries from C. imitator (samples with less than 0.5× coverage were excluded). Among the white-faced capuchin samples, individuals from northern (dry forest) and southern (wet forest) regions form the primary split; secondary splits reflect the individuals from different sites within regions. The short branch lengths of the outgroups are a result of only polymorphic positions within C. imitator being used to construct the tree.

Fig. 6.

(A) Density plot of 1-Mb windows with a slide of 200 kb in northern and southern populations. The distribution of windows from the northern population indicates lower heterozygosity than the southern distribution. The individuals from the southern population show consistently higher values. (B) Long runs of homozygosity in the five largest scaffolds. Blue dots represent windows with depleted heterozygosity. The individuals with the longest runs of homozygosity come from the northern population. (C) PSMC plot of effective population size over time.

Of the 25 capuchin DNA samples (comprising 23 individuals) that we sequenced and mapped to our reference genome, 16 were fecal-derived (). This includes a unique—to our knowledge—high-coverage (12.2×) whole-mammalian genome generated from a fecal sample. When comparing the high-coverage tissue-derived genome from the Santa Rosa site to those generated from our application of fecalFACS, we observed no substantial difference in quality, coverage, heterozygosity, or GC content (). FecalFACS yielded high mapping (median 93%) and low duplication rates (median 6%) (Fig. 2; see details in ). Population clustering relationships were not perturbed by depth of coverage, or source material (tissue-based vs. fecalFACS genomic libraries) (Fig. 5 and ).

Discussion

Local Adaptation to Seasonal Environments.

Patterns of genetic differentiation are consistent with our hypothesis that capuchins inhabiting the tropical dry forests have undergone local adaptation to the extreme seasonality of rainfall and food availability. We observed population-specific variation in genes implicated in water/salt balance, kidney function, muscle wasting, and metabolism. SERPINC1 may reflect selection for drought resistance in the dry-forest capuchins. It has previously been linked to water/salinity balance in mammals; for example, it is overexpressed in the renal cortex of Dahl salt-sensitive rats fed a high salt diet (34). Additionally, the associations of BCAS3 with glomerular filtration rate and AXDND1 with nephrotic syndrome, a kidney disorder that commonly presents with edema and proteinuria, are consistent with this hypothesis. Regarding metabolism, GLIS3 is the only gene other than INS (insulin) known to be associated with all three forms of diabetes (42), and IRF9 has been shown to have a strong effect on body weight in mice, likely through its key role in glucose and lipid metabolism (43). Cautiously, we also point to ABCC5, which nearly missed our FST threshold, but does contain a fixed nonsynonymous SNV in our two populations. The expression of ABCC5 “plays a central role in energy metabolism in mammals,” affecting fat mass and insulin sensitivity in knockout mice and the onset of diabetes in humans (58). Additionally, the significant overabundance of genes associated with diabetic nephropathy in high FST windows is consistent with there being different functional pressures on kidney and metabolic function in these two populations.

That multiple genes in high FST are associated with muscle wasting and serum creatinine and creatinine kinase levels (ITGA7, ISPD, SYNE2, and BCAS3) is of particular interest. Creatinine is a byproduct of the metabolism of creatine phosphate in skeletal muscle, which is normally filtered by the kidneys. It has been used as a clinical biomarker of kidney function, chronic kidney disease (36), and as a monitor of muscle mass (59). Previous field study of urinary creatinine levels of wild capuchins from the northern population of Sector Santa Rosa (SSR) document expected relationships with decreased muscle mass during periods of seasonally low fruit availability, and highlight the biological impacts of seasonal food shortages in the dry forest (16). Study of seasonal variation in SSR capuchin gut microbes (60, 61) also provides evidence of periodic episodes of resource stress and poorer health that differs markedly from those of capuchins inhabiting nonseasonal forests (62). Finally, a recent study has documented mass infant mortality in SSR capuchins following years with longer-than average droughts (15). Together, these observations fit with the notion that animals living in seasonal environments, or pursuing seasonal migrations, are more likely to have weight fluctuations through binge-subsist cycles that map onto food abundance (63). We take this as promising evidence that these genes have been under selection in the northwestern population.

Given that selection operates on both gene function and regulation, we suspect the observed variation is affecting gene expression or enzymatic efficiency, which offers a promising avenue for future research. In addition, overlap between metabolic and immune functions in IRF9 could also suggest inflammation as another target of selection. Because challenging environments will exert selective pressures on numerous other health-related bodily functions, including immunity, future study in this area may be instructive (64). In this vein, further research could examine TLR4, which has been argued to have a key relationship among nutrition, lipids, and inflammation, because of its involvement not only with cytokine expression, but also glucose uptake, body mass, and insulin metabolism (65). Although TLR4 fell outside our high FST window threshold, we did identify a fixed nonsynonymous SNV between populations.

Turning to sensory systems, we found an intact short-wavelength (blue)–sensitive opsin gene OPN1SW, and three alleles of the long-wavelength–sensitive opsin gene (OPN1LW) that are consistent with previous reports for Cebus (66–68). Like most other primates in the Americas, capuchins possess a polymorphic color vision system, characterized by extensive intraspecific genotypic and phenotypic variation (69, 70). This widespread, persistent variation is a fascinating example of balancing selection that has captivated biologists for decades (67, 70, 71). That we find no evidence of differences in the number opsin alleles—nor novel variants—between dry forest and rainforest populations, suggests that these habitat differences are not selecting for variation in color vision intraspecifically. However, we did observe population-specific variation in CCDC66 that could be associated with scotopic (rod-driven) photoreceptor response. Electroretinography of CCDC66−/− mice reveals a significant reduction in scotopic photoreceptor response (44), indicating a potential effect on vision in low-light conditions. Curiously, CCDC66 may impact smell as well as vision. CCDC66−/− mice also display neurodegeneration of the olfactory bulb, and have reduced odor discrimination performance of lemon smells (72). The ecological significance of this result, if any, is unclear at present but may warrant future attention. Turning to chemosensory receptor genes, we find the numbers of functional OR, VR, and TASRs are similar to, or slightly higher than, the number identified in other anthropoid primates (73–76). Interestingly, the VR gene repertoire of capuchins highlights the persistent role of the VR organ that is used in social communication of other mammals and is likely important to the social relationships of capuchins, which frequently engage in urine washing and “hand-sniffing” behaviors (77). The functional significance of the TASR variants we identified is unknown, but may be revealed via cellular expression systems in future research (78). The tuning of the chemosensory receptors to specific stimuli could vary between habitat types and may be elucidated by future work.

A limitation of our study is the relatively small number of sampled individuals from few sampling sites, and that deeper sampling would reduce the chance of spurious results. In order to minimize the possibility of such random hits, we held to strict requirements when identifying evidence of local adaptation. We required that a function directly related to one of our a priori hypotheses be identified as significantly associated by ToppFunn in the strict set of 39 genes that were in the top 0.5% of FST windows and contained a SNV with high or moderate effect (e.g., ITGA7, ISPD, SERPINC1, IRF9, GLIS3, CCDC66) before examining the secondary gene sets (299 in high FST regions and 26 with fixed SNVs) for supplementary support (e.g., BCAS3, ABCC5, TLR4).

Among primates, capuchin monkeys are known for their relatively large brains, derived cognition, and sensorimotor intelligence (2, 5, 6, 9). Accordingly, it is perhaps unsurprising to see positive selection in the Cebus lineage and evidence of shifts in gene function linked to brain function and development relative to other primates. In particular, positive selection in WDR62, BPTF, BBS7, NUP133, and MTOR, and PHF8 is consistent with observations that the capuchin lineage has undergone adaptation linked to brain development. The association of several of these genes with size related brain malformations (46–49), such as microcephaly, suggests that they could be influencing the large relative brain size of capuchins. Furthermore, the evidence of selection in PHF8, which is associated with human cognition, aligns with the link between brain size and intelligence that has been observed in other primates (79). While we highlight here the putative functional roles of these genes, which are based on clinical studies and comparative genomics, we acknowledge that further examination of their function in the context of capuchin biology is warranted.

In the context of longevity, it is noteworthy that we observed genes under selection associated with DNA damage response, metabolism, cell cycle, and insulin signaling (13). Of particular interest are PARP1, MTOR, SREBF1, INSR1, and HTT. Damage to the DNA is thought to be a major contributor to aging (80). Previous studies have also shown that genes involved in DNA-damage responses exhibit longevity-specific selection patterns in mammals (14). It is therefore intriguing that PARP1, a gene suggested to be a determinant of mammalian aging (81), is under selection in capuchins. Another large body of research has associated the MTOR with aging and longevity in various organisms (82), making it a prime candidate for therapeutic interventions in aging (83). Other genomes of long-lived mammals also revealed genes related to DNA-repair and DNA-damage responses under selection (84, 85). Intriguingly, short-lived species also exhibit genes under selection related to insulin receptors, raising the possibility that the same pathways associated with aging in model organisms are involved in the evolution of both short and long lifespans (86), an idea supported by our results. Of course, because aging-related genes often play multiple roles, for example in growth and development, it is impossible to be sure whether selection in these genes is related to aging or to other life-history traits, like growth rates and developmental times, that in turn correlate with longevity (87). Therefore, although we should be cautious about the biological significance of our findings, it is tempting to speculate that, like in other species, changes to specific aging-related genes or pathways, could contribute to the longevity of capuchins. Additionally, it is difficult to identify the evolutionary time depth beyond the Saimiri/Cebus branch [estimated time: 16.07 Mya; confidence interval (CI): 15.16 to 17.43 Mya (88)] based on the present study. Future analyses of other capuchin species, in both extant genera (Cebus and Sapajus) along with other Cebids and platyrrhines will better reveal more recent selection. Future research assessing structural variation across the genome could also be fruitful, particularly in cases such as the copy number of CD33 related Siglec genes, which have previously been linked to longevity in mammals (89).

Finally, it is worth noting that subtle changes in the duration of expression of regulatory genes could have profound effects on both lifespan and brain development. Evidence of the temporal importance of the onset of neurogenesis and effects on interspecific variation in neuron composition are present in comparative study of owl and capuchin monkeys, for example (90). These changes are not possible to detect with the orthologous gene positive-selection scan approach that we take here; future research examining likely differences in neuroembryological parameters, along with variation in gene-expression profiles across the lifespan, would complement and extend these contributions.

Population Genomics with FecalFACS.

We observe a clear demarcation in population structure between the northern dry-forest and southern wet-forest populations in Costa Rica. Higher levels of heterozygosity in the south and lower levels in the northwest are in accordance with the current hypothesis that capuchins dispersed northwards across Costa Rica. White-faced capuchins are the most northerly distributed member of the Cebinae, having dispersed over the isthmus of Panama in a presumed speciation event with Cebus capucinus in South America ∼1.6 Ma (91, 92). After expanding during the early late Pleistocene, white-faced capuchins appear to have undergone a dramatic reduction in effective population size. This pattern predates the movement of humans into Central America, and could reflect a series of population collapses and expansions caused by glacial shifts and fluctuating forest cover availability during the Pleistocene. C. imitator in SSR are near the northernmost limits of their range, which extends as far north as Honduras, and may represent some of the least genetically diverse members of their species. Given the limitations of the available sampling sites, it is possible that the appearance of an ecological divide is actually evidence of isolation by distance; however, given that the single individual from Cañas clusters closely with the individuals from SSR, despite a geographic distance of more than 100 km, we suggest that isolation by distance does not completely explain the population differentiation.

Through a use of flow cytometry/FACS, we successfully generated a high-coverage, minimally biased mammalian genome solely from feces, along with a low-coverage SNP dataset that is suitable for population assignment and clustering. We also used the high-coverage fecalFACS genome to call SNPs (along with traditional blood and tissue genomes), demonstrating that noninvasively collected fecal samples can successfully be used to generate primate genomic resources. The clustering patterns in our trees and PCA plots do not reveal any samples that deviate from their expected geographic or ecological origin. These relationships are robust to both the coverage levels (<1× to >50×) and biological origins (feces, tissue, and blood) of the samples. The tight geographic clustering of individuals within the SSR sampling locale provides reasonable evidence that there is no substantial effect from fecalFACS on population structure. One possible limitation of fecalFACS could be the coisolation of cells from prey items in carnivorous/omnivorous mammals. On rare occasions, capuchins eat small mammalian prey, such as infant coatimundis and squirrels (93), and we did observe a small contribution (<2% of reads in two samples) that could be from mammalian prey (). However, we successfully removed these reads with BBsplit (94) mapping (), and we did not notice any effect of contamination in our results (Fig. 5). We suspect that any dietary contribution is low, as partially digested cells would be more likely to be degraded in a way that further distinguishes them in size and granularity (i.e., the forward- and side-scatter parameters used in our FACS gating) from epithelial cells sloughed from the lining of the colon. Given that fecalFACS is cost-effective and minimizes the biases that commonly occur in traditional bait-and-capture approaches to the enrichment of endogenous DNA from feces, and does not require costly, impractical preservation of biomaterial in liquid nitrogen but rather uses room-temperature stable storage in RNAlater, this method offers great benefits to the field of mammalian conservation and population genomics. Our approach is promising for future studies of wild animals where invasive sampling is not feasible or ethical.

Summary

In our population level analysis of wet- and dry-forest capuchins, we observed both evidence of population structure between and local adaptation to these different habitats. In particular, we identified selection in genes related to food and water scarcity, as well as muscular wasting, all of which have been observed during seasonal extremes in the dry forest population. Interspecific analyses revealed evidence of selection in C. imitator on genes involved with brain development and longevity. These results are in accordance with the remarkably long life span, large brain, and high degree of sensorimotor intelligence that has been observed in capuchins. These genes are good candidates for further investigation of traits that have evolved in parallel in apes and other mammals. Importantly, our insights were made possible through an annotated reference assembly of a capuchin monkey and a use of flow cytometry/FACS to isolate epithelial cells from mammalian feces for population genomics. FecalFACS allowed us to generate both a high-coverage, minimally biased mammalian genome solely from feces, as well as low-coverage SNV datasets for population level analyses.

Methods

Study Populations and Sample Collection.

Central American white-faced capuchins (Cebus imitator), a member of the gracile radiation of capuchins (genus Cebus), were recently recognized as a species, distinct from C. capucinus in South America (92). C. imitator occupies a wide diversity of habitats, spanning lowland rainforests and cloud forests in Panama and throughout southern, eastern, and central Costa Rica, and tropical dry forests in northwestern Costa Rica and Nicaragua. The annual precipitation and elevation of rainforest versus dry forest biomes in their current ranges vary dramatically, leading to considerable variation in the resident flora and fauna (95, 96). We sampled individual Costa Rican capuchin monkeys from populations inhabiting two distinct habitats: 1) Lowland rainforest around Quepos, Puntarenas Province; and 2) tropical dry forest at two sites in Guanacaste Province. In total, we collected samples from 23 capuchins, a list of which is provided in .

We sampled capuchins inhabiting a lowland tropical rainforest biome by collaborating with Kids Saving the Rainforest (KSTR) in Quepos, Costa Rica. We acquired blood samples from four wild capuchins from nearby populations who were undergoing treatment at the facility (although we were unable to collect paired fecal samples). For one of these individuals, an adult male white-faced capuchin that was mortally wounded by a vehicle in Costa Rica, we additionally sampled tissues from several organs. DNA derived from the kidney was used for the reference genome assembly.

We collected 21 samples from 19 individuals in the northern tropical dry forest. Sixteen fecal samples and 4 tissue samples were from free-ranging white-faced capuchin monkeys (C. imitator) in SSR, part of the Área de Conservación Guanacaste in northwestern Costa Rica, which is a 163,000-ha tropical dry forest nature reserve (Fig. 1). Behavioral research of free-ranging white-faced capuchins has been ongoing at SSR since the 1980s, which allows for the reliable identification of known individuals from facial features and bodily scars (97). The 16 fresh fecal samples were collected from 14 white-faced capuchin monkeys immediately following defecation (). We placed 1 mL of feces into conical 15-mL tubes prefilled with 5 mL of RNAlater. RNAlater-preserved fecal samples were sent to the University of Calgary, where they were stored at room temperature for up to 3 y. To evaluate other preservation methods, we also collected two additional capuchin monkey fecal samples (SSR-FL and a section of SSR-ML), which we stored in 1× PBS buffer and then froze in liquid nitrogen with a betaine cryopreservative (98). Given the logistical challenges of carrying liquid nitrogen to remote field sites, we prioritized evaluation of samples stored in RNAlater. We also collected tissue and blood samples opportunistically. During the course of our study, four individual capuchin monkeys died of natural causes at SSR, from whom we were able to collect tissue samples, which were stored in RNAlater. Additionally, we collected a blood sample from one northern dry-forest individual housed at KSTR that originated near the town of Cañas, which is ∼100 km southwest of SSR.

Summary of Genome-Wide Sequencing, Genome Assembly, and Gene Annotation.

We assembled a reference genome for C. imitator from DNA extracted from the kidney of a male Costa Rican individual (KSTR64) using a short-read approach (Illumina HiSeq 2500) with ALLPATHS-LG. To improve the quality of gene annotation, we isolated total RNA from the whole blood of an adult male white-faced capuchin (ID: CNS-HE) permanently residing at the KSTR wildlife rehabilitation center. The capuchin genome assembly was annotated with the National Center for Biotechnology Information (NCBI) pipeline previously described in https://www.ncbi.nlm.nih.gov/books/NBK169439/.

Our reference genome assembly for C. imitator is composed of 7,742 scaffolds (including single contig scaffolds) with an N50 scaffold length of 5.2 Mb and an N50 contig length of 41 kb. The final ungapped assembly length is 2.6 Gb (GenBank accession: GCA_001604975.1). Our estimate of total interspersed repeats using WindowMasker (99) output is 45.8%. The numbers of annotated genes are 20,740 and 9,556 for protein-coding and noncoding genes, respectively (). Measures of gene representation using the known human RefSeq set of 56,230 transcripts show an average of >94% coverage with a mean identity of 92.5%. Overall, our draft assembly metrics and gene representation are consistent with other nonhuman primate short-read reference assemblies (100). See for further methodological detail.

Summary of Positive Natural Selection Analysis through Codon-Based Models of Evolution and Enrichment Tests.

The phylogenetic arrangement in this study included 14 species as outgroups to C. imitator: Three Platyrrhini (Callithrix jacchus, Aotus nancymaae, Saimiri boliviensis), six Catarrhini (Macaca mulatta, Rhinopithecus roxellana, Nomascus leucogenys, Pan troglodytes, Homo sapiens, Pongo abelii), one Strepsirrhini (Microcebus murinus), one rodent (Mus musculus), and three Laurasiatheria (Canis lupus familiaris, Bos taurus, and Sus scrofa). We identified 7,519 Orthologous Groups (OGs) present among the 15 species. We identified 612 genes under positive selection (P < 0.05 after false-discovery rate correction) in the Cebus lineage using the branch model of codeml in PAML. We also performed a branch-site test and identified a second set of 746 genes under positive selection in Cebus (Dataset S1). We performed functional annotation analysis using ToppFunn in the ToppGene Suite (45) with default parameters. To ascertain which ontology processes the genes with signals of positive selection were involved, we focused the enrichment analysis on two functional categories: GO Biological Processes (BP) and the DisGeNET BeFree disease database. Finally, the genes with positive selection signal were intersected with the GenAge and CellAge databases (build 19,307 genes) (53). See for further methodological detail.

FecalFACS.

Before isolating cells by FACS, fecal samples were prepared using a series of washes and filtration steps. Fecal samples were vortexed for 30 s and centrifuged for 30 s at 2,500 × g. Then the supernatant was passed through a 70-μm filter into a 50-mL tube and washed with Dulbecco's phosphate-buffered saline (DPBS). After transferring the resultant filtrate to a 15-mL tube, it was centrifuged at 1,500 rpm for 5 min to pellet the cells. Then we twice washed the cells with 13 mL of DPBS. We added 500 μL of DPBS to the pellet and refiltered through a 35-μm filter into a 5-mL FACS tube. We prepared a negative control (to control for auto-fluorescence) with 500 μL of DPBS and one drop of the cell solution. To the remaining solution, we added 1 μL of AE1/AE3 Pan Cytokeratin Alexa Fluor 488 antibody (ThermoFisher: 53-9003-82) or TOTO-3 Iodide (642/660) DNA stain, (ThermoFisher T3604), which we allowed to incubate at 4 °C for at least 30 min.

We isolated cells using a BD FACSAria Fusion (BD Biosciences) flow cytometer at the University of Calgary Flow Cytometry Core. To sterilize the cytometer’s fluidics before processing each sample, we ran a 3% bleach solution through the system for 4 min at maximum pressure. We assessed background fluorescence and cellular integrity by processing the negative control sample prior to all prepared fecal samples. For each sample we first gated our target population by forward- and side-scatter characteristics that were likely to minimize bacteria and cellular debris (). Secondary and tertiary gates were implemented to remove cellular agglomerations. Finally, we selected cells with antibody or DNA fluorescence greater than background levels. In cases when staining was not effective, we sorted solely on the first three gates. Cells were pelleted and frozen at −20 °C. A detailed fecalFACS laboratory protocol is available in .

We isolated a median of 2,206 cells, with a range of 129 to 62,201 (). The total amount of DNA per sample was low, ranging from 2.96 to 21.50 ng, with a median value of 7.85 ng, but in each case allowed for the successful generation of a sequencing library (). The number of cells was not significantly correlated with the amount of extracted DNA (R = 0.227; 95% CI [−0.345, 0.676]; t = 0.808, P > 0.05) or mapping rate (R = −0.204; 95% CI [−0.663, 0.367]; t = −0.721; P > 0.05). Median mapping rates reached 93% (range: 55 to 98%) with BWA-MEM and 82% (range: 11 to 95%) with the more stringent BBsplit settings (Fig. 2 and ). Read duplication levels were low, with a median value of 6% (range: 2 to 40%) resulting in 63% (range: 8 to 92%) of reads being unique and mapping to the C. imitator 1.0 genome. The amount of duplicate reads was distributed bimodally across individuals, with reads from five samples having substantially higher duplication rates than the remaining nine. The rate of duplication was significantly correlated (R = −0.751; 95% CI [−0.917, −0.366]; t = −3.94; P < 0.01) with the number of cells (log10-transformed), decreasing sharply above a threshold of about 1,000 cells (Fig. 2). By sorting fecal samples with FACS, we substantially increased the percentage of reads mapping to the target genome. We selected five samples at random (SSR-CH, SSR-NM, SSR-LE, SSR-PR, SSR-SN) to compare pre- and post-FACS mapping rates. The mapping rates of unsorted feces ranged from 10 to 42%, with a median of 14% (Fig. 2). After flow sorting aliquots of these fecal samples, we obtained significantly higher mapping rates (V = 15, P < 0.05) for each sample, ranging from 64 to 95%, with a median of 85%, resulting in a median 6.07-fold enrichment.

DNA Extraction and Shotgun Sequencing.

We extracted fecal DNA (fDNA) with the Qiagen DNA Micro kit, following the “small volumes of blood” protocol. To improve DNA yield, we increased the lysis time to 3 h, and incubated 50 µL of 56 °C elution buffer on the spin column membrane for 10 min. DNA concentration was measured with a Qubit fluorometer. Additionally, to calculate endogenous DNA enrichment, we extracted DNA directly from five fecal samples prior to their having undergone FACS. We extracted DNA from the nine tissue and blood samples using the Qiagen Gentra Puregene Tissue kit and DNeasy blood and tissue kit, respectively.

For the fecal samples, DNA was fragmented to 350 bp with a Covaris sonicator. We built whole-genome sequencing libraries with the NEB Next Ultra 2 kit using 10 to 11 PCR cycles. Fecal genomic libraries were sequenced on an Illumina NextSeq (2 × 150 PE) at the University of Calgary genome sequencing core and an Illumina HiSeq 4000 at the McDonnell Genome Institute at Washington University in St. Louis. Using half of one HiSeq 4000 lane, we achieved an average coverage of 12.2× across the C. imitator 1.0 genome (sample SSR-ML). Other fecal samples were sequenced to average depths of 0.1 to 4.4× (). High-coverage (10.3 to 47.6×), whole-genome shotgun libraries were prepared for the blood and tissue DNA samples and sequenced on an Illumina X Ten system at the McDonnell Genome Institute. For population analyses within capuchins, we mapped genomic data from all 23 individuals sequenced () to the reference genome.

Mapping and SNV Generation.

We called SNVs for each sample independently using the Cebus imitator 1.0 genome and the GATK UnifiedGenotyper pipeline (see for further methodological detail). We included reads from all nine tissue/blood samples and one frozen fecal sample with high coverage (SSR-ML). In total, we identified 4,184,363 SNVs for downstream analyses. To remove potential human contamination from sequenced libraries, we mapped trimmed reads to the C. imitator 1.0 and human (hg38) genomes simultaneously with BBsplit. Using default BBsplit parameters, we binned separately reads that mapped unambiguously to either genome. Ambiguously mapping reads (i.e., those mapping equally well to both genomes) were assigned to both genomic bins, and unmapped reads were assigned to a third bin. We calculated the amount of human genomic contamination as the percentage of total reads unambiguously mapping to the human genome (). After removing contaminant reads, all libraries with at least 0.5× genomic coverage were used for population structure analysis.

In order to test the effect of fecalFACS on mapping rates, we selected five samples at random (SSR-CH, SSR-NM, SSR-LE, SSR-PR, SSR-SN) to compare pre- and post-FACS mapping rates. To test for an increase in mapping percentage, we ran a one-sample paired Wilcoxon signed-rank test on the percentages of reads that mapped exclusively to the Cebus genome before and after FACS. Additionally, we ran Pearson’s product moment correlations to test for an effect of the number of cells (log10-transformed) on rates of mapping, read duplication, and nanograms of input DNA. The above tests were all performed in R.

High-Coverage Fecal Genome Comparison.

We made several comparisons between our high-coverage feces-derived genome and the blood/tissue-derived genomes using window-based approaches. For each test, the feces-derived genome should fall within the range of variation for members of its population of origin (SSR). Deviations from this, for example all fecal genomes clustering together, would indicate biases in our DNA isolation methods. To assess this, we constructed 10-kb windows with a 4-kb slide along the largest scaffold (21,314,911 bp) in the C. imitator reference genome. From these windows, we constructed plots of coverage density and the distribution of window coverage along the scaffold. Secondly, we assessed the level of heterozygosity in 1-Mb/200-kb sliding windows throughout the ten largest scaffolds. For each high-coverage genome, we plotted the density distribution of window heterozygosity. We measured genome-wide GC content with the Picard Tools CollectGcBiasMetrics function. The percentage of GC content was assessed against the distribution of normalized coverage and the number of reads in 100-bp windows per the number reads aligned to the windows.

Population Structure.

Given the large degree of difference in coverage among our samples, (less than 1× to greater than 50×), we performed pseudodiploid allele calling on all samples. For each library, at each position in the SNV set, we selected a single, random read from the sequenced library. From that read, we called the variant information at the respective SNV site for the given library. In doing so, we generated a variant call format (VCF) file with a representative degree of variation and error for all samples.

To assess population structure and infer splits between northern and southern groups of Costa Rican white-faced capuchins, we constructed principal components plots with EIGENSTRAT (101) and built population trees with TreeMix (102). Because we ascertained variants predominantly with libraries that were of tissue/blood origin, we built PCs solely with SNVs from these libraries and projected the remaining fecal libraries onto the PCs. For our maximum-likelihood trees, we used two outgroups (Saimiri sciureus and Cebus albifrons), with S. sciureus serving as the root of the tree. Given the geographic distance and anthropogenic deforestation between northern and southern populations, we assumed no migration. To account for linkage disequilibrium, we grouped SNVs into windows of 1,000 SNVs. Population size history was inferred from the highest coverage nonreference individual, SSR-RM08 with PSMC, using default parameters (https://github.com/lh3/psmc).

Local Adaptation, FST, Heterozygosity, Relatedness.

For all analyses of local adaptation and heterozygosity between populations, we excluded individuals from our low-coverage dataset. We tested for the degree of relatedness among all high- and low-coverage individuals using READ (103) and identified two of the high-coverage individuals from SSR, SSR-ML and SSR-CR, as potential first-degree relatives (). For all statistical analyses of high-coverage samples, we removed SSR-ML, because SSR-CR was sequenced to higher average depth.

For each individual, we calculated heterozygosity in 1-Mb/200-kb sliding windows across the genome for all scaffolds at least 1 Mb in length. Windows were generated with BedTools windowMaker (104) and heterozygosity was calculated as the per site average within each window. Based upon a visual inspection of the average heterozygosity values across the genome (Fig. 6), we classified a window as part of a run of homozygosity if the window’s average heterozygosity fell below 0.0002. Descriptive statistics and two-sided Wilcoxon tests were calculated in R.

For each high-coverage sample, we calculated the Hudson’s FST ratio of averages (105) in 20-kb windows with a slide of 4-kb across the genome. Among the genes present in each window in the top 0.5% and top 0.1% of FST values, we searched for SNPs with high or moderate effects using SnpEff and identified those SNPs with high FST values (> 0.75) using VCFtools. We searched for functional enrichment of our population gene set using ToppFun in the ToppGene Suite (45). In an effort to identify candidate genes for further investigation, we set low thresholds for significance (false-discovery rate P < 0.1 and the minimum number of genes per category to 1).

Chemosensory Genes.

The chemosensory behaviors of capuchins have been well-studied (73, 75), and taste and olfaction are suspected to play important roles in their foraging ecology; 614 orthologous olfactory receptor genes were identified in the Cebus reference genome using previously described methods (106). Briefly, putative OR sequences were identified by conducting TBLASTN searches against the capuchin reference assembly using functional human OR protein sequences as queries with an e-value threshold of 1e-20. For each reference-derived OR gene, we added 500 bp of flanking sequence to both the 5′ and 3′ ends in a bed file. For each individual, we extracted the OR gene region from the gVCF and generated a consensus sequence defaulting to the reference allele at variable site using bcftools (107). The number of intact, truncated, and pseudogenized OR genes in each individual were identified using the ORA pipeline (108). We considered an OR gene to be putatively functional if its amino acid sequence was at least 300 amino acids in length. ORA was further used to classify each OR into the appropriate class and subfamily with an e-value cutoff of 1e-10 to identify the functional OR subgenome for each individual (108). TASR and VR genes were identified in the NCBI genome annotation, and variable positions were located by scanning the VCF with VCFtools and bash. Cebus opsin tuning site positions have been identified previously (68). With the high-coverage dataset, we identified the alleles present at each locus in the VCF. For each low-coverage fecal-derived genome, we located the position of the tuning site in the bam file using SAMtools tview and manually called the variant when possible.

Ethics Approval and Consent to Participate.

This research adhered to the laws of Costa Rica, the United States, and Canada and complied with protocols approved by the Área de Conservación Guanacaste and by the Canada Research Council for Animal Care through the University of Calgary’s Life and Environmental Care Committee (ACC protocol AC15-0161). Samples were collected with permission from the Área de Conservación Guanacaste (ACG-PI-033-2016) and CONAGEBIO (R-025-2014-OT-CONAGEBIO; R-002-2020-OT-CONAGEBIO). Samples were exported from Costa Rica under permits from CITES and Área de Conservación Guanacaste (2016-CR2392/SJ #S 2477, 2016-CR2393/SJ #S 2477, DGVS-030-2016-ACG-PI-002-2016; 012706) and imported with permission from the Canadian Food and Inspection Agency (A-2016-03992-4).