These results, although preliminary, are in line with the evidence from other studies indicating that even opportunistic dispersal can be sufficient to prevent inbreeding, as long as it is unconstrained by habitat fragmentation or other factors. The importance of unconstrained dispersal for inbreeding avoidance was supported in a simulation study based on empirical dataset from golden-crowned sifakas, Propithecus tattersalli, showing that high levels of outbreeding can be maintained in a population by a combination of social structure and unconstrained dispersal but without the need for active inbreeding avoidance mechanisms 68. The link between dispersal and inbreeding risk was further indicated by studies demonstrating a correlation between dispersal distances and inbreeding level (e.g., in great tits, Parus major: Szulkin and Sheldon, 2008). At our study site, the habitat was undisturbed, and dispersal was most likely unconstrained, ensuring passive inbreeding avoidance. As indicated by one case where pair mates were second-degree kin, occasional inbreeding can still occur in such populations and is presumably tolerated.
In addition to dispersal and preferential mating with unrelated individuals, another potential way to avoid inbreeding is through engaging in EPC. Positive relationship between EPP rates and pair mate relatedness was demonstrated in many bird species 7,10,11, but pair-living mammals do not seem to use this strategy often. In pair-living mammals, a similar strategy was, to our knowledge, only demonstrated in one species, fat-tailed dwarf lemur, where females sharing more MHC-supertypes with their social partner were shown to engaged in more EPC 14. In our study population, the absence of evidence for EPP further confirms our suggestion that dispersal in this habitat is unconstrained and the potential for inbreeding is low, rendering EPC unnecessary.
Summing up, the current study is the first to examine the link between mating system, mate choice and dispersal in a wild population of a pair-living primate. We showed that coppery titis in our study population are mostly genetically monogamous, likely as a result of a strong pair bond enabling effective mate guarding and relatively low population density limiting the opportunities for extra-pair copulations. We further showed that coppery titis, despite exhibiting no active inbreeding avoidance via mate choice, still had low relatedness between pair mates. Our results indicate that even opportunistic dispersal, as long as it is not constrained, can create sufficient genetic dissimilarity between opposite sexes to render active mate choice and extra-pair copulations unnecessary. Future studies with larger sample size will be needed to examine the extent of genetic monogamy in different populations of coppery titis and to further investigate dispersal patters. In particular, to examine if titis indeed lack the mechanisms for active inbreeding avoidance via kin discrimination, it will be necessary to compare mating patterns and levels of inbreeding in continuous vs. fragmented or isolated populations. Finally, the absence of relatedness- and heterozygosity-based mate choice in our study population, of course, does not mean that mate choice does not occur in titis at any level. To better understand mating patterns in titis, future studies will have to examine if mate choice is based on other factors, such as, e.g., variation in MHC loci, body condition or the size or quality of the territory.
Methods
Study site and study population
The study was conducted at the Estación Biológica Quebrada Blanco in the north-eastern Peruvian Amazon (4°21’S, 73°09’W) in June 2017 – September 2019. Study individuals belonged to 14 family groups (Supplementary Table S1, Fig. 1), seven of which (Groups 1–7) were also subject to behavioral studies conducted in June–December 2017 and June–December 2018  38,70,71. Between the periods of behavioral data collection, the groups were monitored for 2–3 days per month, and genetic samples were collected continuously from the beginning of the study until September 2019. Group 1 had been habituated to the presence of human observers and studied intermittently since 1997; the other groups were habituated during this study. We individually identified all the study animals based on the combination of body size, tail shape and colouration, genitalia shape, and natural marks.
Titis typically give birth to a single infant once a year 33,57,72. In our study population, most of the births occurred between October and February and only one occurred in June (Supplementary Table S1). As the offspring disperse at the age of 2–4 years 33,57,72, the pedigree in our study comprised up to 5 generations of offspring per group (Supplementary Table S1).
Fecal sample collection and DNA extraction
We collected fecal samples from 41 individuals (3-15 samples per individual) living in 14 family groups, including 18 putative offspring of 9 family groups (1 to 5 offspring per group). Five other groups either did not have offspring during the study period (or they had disappeared before we could collect samples) or the samples could not be collected because the offspring were still very young and thus their defecations too diminutive to be detected. Also due to differential habituation to the presence of humans, for some groups we could not obtain samples from all group members. For those groups that were habituated in the beginning of the study period, we collected samples from offspring from several consecutive years.
Fecal samples were collected immediately after an identified individual was seen defecating. We dried the samples in 15 mL falcon tubes containing silica gel beads (Carl Roth, Karlsruhe, Germany) and stored them at ambient temperature, replacing the silica beads when necessary, until shipping to Germany.
We extracted DNA (at least two samples per individual for all animals except one offspring of Group 10; see Results for more details) from ca. 200 mg of feces using the Macherey-Nagel NucleSpinã DNA stool kit with a final elution of the DNA in 50 mL elution buffer. DNA concentration of the extracts was measured using a NanoDrop Spectrophotometer (ND-1000, PEQLAB Biotechnologie GmbH) and a Qubit Fluorometer (Thermo Fisher).
Microsatellite genotyping
As published microsatellite loci for titi monkeys 73–75 revealed unreliable results for our study species, we established a new set of 27 di-repeat microsatellite loci that can be universally applied to New World monkeys (details are described in Supplementary Materials and Supplementary Tables S3–4). To simplify library preparation for genotyping by sequencing, we added adapter nucleotide sequences to the 5’ end of the locus-specific primers.
We amplified all 27 loci in a single multiplex PCR using the Qiagen Multiplex PCR Kit (Qiagen) with a total volume of 25 mL and containing 12.5 mL 2x Multiplex Master Mix, 1 mL of primer pool (0.2 mM of each primer), 1 mL of DNA extract (ca. 200 ng total DNA) and 10.5 mL of RNase-free water. Amplifications were performed with initial denaturation at 95°C for 15 min, 40 cycles of denaturation at 94°C for 30 sec, annealing at 57 °C for 1.5 min, extension at 72 °C for 1 min and a final extension at 60°C for 30 min. PCR products were checked on 1.5% agarose gels together with non-template controls. To prevent false homozygosity due to allelic dropout, we repeated each multiplex reaction three times 76. In some samples, the total multiplex reaction with all 27 loci yielded low number of sequencing reads; in these cases, we additionally amplified the loci in three separate multiplex reactions with the following primer pools: chr01b­–chr07a, chr08a–chr12a, chr12b­–chrXa, as this method usually yielded more reads (see Supplementary Materials for details). The reactions and PCR conditions for three separate multiplex reactions were the same as for the total multiplex reaction.
Following amplification, we pooled 5 mL of each multiplex PCR product (or of each PCR product of three separate multiplex reactions), purified the pooled products with the Monarch PCR & DNA Cleanup Kit (New England BioLabs) and quantified them using Qubit Fluorometer (Thermo Fisher). To prepare sequencing libraries, we performed indexing PCRs using Kapa HiFi Hotstart ReadyMix PCR Kit (Roche) with a total volume of 25 mL containing 12.5 mL 2x Kapa HiFi Hotstart ReadyMix, 1 mL (0.5 mM) of each indexing primer (containing individual barcodes) and 100 ng purified PCR product. Indexing PCRs were done with an initial denaturation step at 98°C for 45 sec, followed by 4 cycles of denaturation at 98°C for 15 sec, annealing at 62°C for 30 sec and extension at 72 °C for 30 sec, and a final extension step at 72°C for 1 min. Full-length libraries were purified with the Monarch PCR & DNA Cleanup Kit (New England BioLabs) and quantified using Qubit Fluorometer (Thermo Fisher). Fragment sizes and molarities were quantified using a Bioanalyzer 2100 (Agilent Technologies). Libraries were diluted to 10 nM and then pooled and sequenced using Miseq Reagent Kit v2 with PhiX DNA (Illumina) added on the MiSeq system (Illumina). Sequencing was performed with 51 forward and 251 reverse read cycles. Only the reverse reads were used for further analysis, while forward reads were only used for MiSeq quality control.
To control for possible misidentification of animals in the field, we genotyped most individuals from 2–3 independent samples. We also used a PCR-based sexing assay 77 to confirm reported sex (and to sex young individuals for whom sex could not be identified in the field). To control for laboratory mistakes, we genotyped each sample twice, leading to at least four genotypes per individual.
After sequencing, the samples were demultiplexed using MiSeq Reporter software and then processed using the CHIIMP analysis pipeline 76. The CHIIMP pipeline calls alleles by first producing unique sequences with relevant attributes (read counts, sequence length, etc.) for each MiSeq sequence file, querying the sequences for potential PCR artifacts, such as stutter sequences, and then removing all sequences that do not match the locus attributes. All alleles called by CHIIMP were manually checked to validate the results and to correct automated allele calling for those loci that contain “wobble” positions in the primer sequences and are incorrectly processed by CHIIMP. We used a cutoff of 250 reads. Additionally, we accepted alleles if they yielded >100 reads in at least three genotypes obtained per individual. Alleles with <100 reads were not called.
Of 27 loci, nine either consistently failed to amplify in our study animals (chr06b, chr11f, chr16b) or proved to be monomorphic (chr02a, chr02b, chr04a, chr10b, chr12a, chr13b) and were excluded from further analysis. The final set consisted of 18 loci, including 17 autosomal and one X-linked locus (chrXa) (Supplementary Table S5). All animals were genotyped at a minimum of 14 loci (16.8 loci on average), and the mean number of alleles per locus was 8.9.
We checked all loci for the presence of null alleles, allelic dropout and stuttering using Micro-Checker 2.2.5 78. We assessed Hardy–Weinberg equilibrium (HWE) and calculated observed and expected heterozygosity with PopGenReport 2.2.2 79. Since the presence of family structure can cause deviations from HWE and bias population genetic analyses, especially in monogamous species, we only included adults in this analysis. The analysis indicated that the population was in HWE. Two loci, chr01b and chr21a, departed from HWE, likely due to the presence of relatives in a study group and/or small sample size.
One of these two loci, chr01b, also showed evidence of null alleles. As the locus did not show any mismatches for the known mother/offspring dyads (see below), we ran all further analyses using two sets of data, one with the full set of loci and another one with locus chr01b excluded. Since the results from these two sets did not differ substantially, we present all further results only for the reduced data set.
Mitochondrial DNA (mtDNA) genotyping
We genotyped all individuals at the hypervariable region I of the mitochondrial control region using primers 5’-TACCTCGGTCTTGTAAACCG-3’ and 5’-AGGTAGGAACCAGATGCCG-3’, newly designed on the basis of mitochondrial genomes of New World monkeys available in GenBank. PCR reactions with a total volume of 30 µl contained 1 U BiothermTaq 5000 (Genecraft), 1x reaction buffer, 0.16 mM of each dNTP, 0.33 µM of each primer and ca. 100ng total DNA. PCRs were performed with initial denaturation at 95 °C for 2 min, 40 cycles of denaturation at 95 °C for 1 min, annealing at 58 °C for 1 min, extension at 72 °C for 1 min and a final extension at 72 °C for 5 min. PCR products were run on 1% agarose gels, excised from the gel and then purified with the Monarch DNA Gel Extraction Kit (New England BioLabs) and sequenced on an ABI 3130xL sequencer using both amplification primers and the BigDye Cycle Sequencing Kit (Applied Biosystems). Sequence electropherograms were checked with 4Peaks 1.8 (https://nucleobytes.com/4peaks/index.html) and manually edited and assembled in SeaView 4.5.4 80; all haplotypes were 567 bp long.
Statistical analyses
As X-linked loci are haploid in males and cannot be treated in the analyses in the same way as autosomal loci, all the following statistical tests were performed using the set of 16 autosomal loci for both sexes. The data for the X-linked locus chrXa was used separately to manually check for allelic mismatches between candidate parents and offspring in the parentage analyses.
(1) Parentage analyses
Parentage was assigned using Cervus 3.0 81. Cervus compares likelihood ratios of the two most likely candidate parents and assigns parentage based on statistical thresholds generated during the simulation analysis. For Cervus analysis, we used a simulation of 100,000 offspring, an error rate of 0.01, 90% relaxed and 95% strict confidence level, and accounted for relatedness of candidate mothers and fathers. Relatedness was calculated with the R package related 1.0 82 using Wang’s estimator r 83. This estimator was chosen because it performed best in simulations, showing the highest correlation between observed and expected values for our set of loci. Additionally, we used Colony 2.0.6.5 84 to verify parentage assignments from Cervus. Unlike Cervus, Colony reconstructs a full pedigree, inferring sibship and parentage among individuals by comparing the likelihood of different clusters of individuals and maximizing group rather than pairwise likelihoods. For this analysis, we used an error rate of 0.01, male and female polygyny, and a sibship size prior of 1.6, calculated as the average number of offspring per family group in our study population.
For both Cervus and Colony analyses, the set of candidate fathers included all sampled adult males plus the oldest subadult male from Group 6 that had dispersed from his natal group in the beginning of the study and could have sired offspring by the end of the sampling period. The set of candidate mothers included all adult females that shared their mtDNA haplotype with candidate offspring. For seven offspring (Supplementary Table S1), the mothers were known because they were seen nursing them. To test the reliability of our parentage estimates, we ran the analyses twice, with and without the respective set of known mother-offspring pairs. Combined non-exclusion probability for the set of 16 autosomal loci (with chr01b excluded) was 9.9x10-5 for the first parent, 3.4x10-7 for the second parent, and 9.0x10-12 for the parent pair, calculated using Cervus.
(2) Relatedness-based mate choice
To test if titis avoid mating with related individuals, we compared relatedness between real and randomly generated mating partners using the pairwise relatedness estimator implemented in STORM 85. First, STORM calculates the relatedness of real mating pairs using the estimator of Li et al. (1993), with each locus weighted using the method described in Lynch and Ritland (1999) and Van de Casteele et al. (2001). Then, the program calculates the expected relatedness of mating pairs if the mating is random with respect to relatedness; this is done by generating mating pairs from female and male breeding pools over 1000 iterations and averaging the relatedness values for each iteration. The obtained distribution is then compared to the averaged relatedness of real mating pairs. Our sample included ten real mating pairs, and the breeding pool consisted of 12 females and 12 males. This included all sampled adults and the oldest subadult male from Group 6.
(3) Heterozygosity-based mate choice
To test if titis show any heterozygosity-based mating pattern, we compared individual heterozygosity levels between pair mates. To estimate individual heterozygosity, we calculated homozygosity by loci (HL), a microsatellite-derived measure that weights the contribution of each locus to the homozygosity value depending on their allelic variability, implemented in R function GENHET 3.1 89. To test if HL is correlated between pair mates, we used a two-tailed Pearson correlation analysis.
(4) Dispersal patterns
To examine if dispersal distances differ between sexes, we compared the diversity of mtDNA haplotypes, relatedness, and heterozygosity among adult females and males. MtDNA haplotype and nucleotide diversity was calculated and compared using a permutation test implemented in R function genetic_diversity_diff 1.0.6 (Alexander et al., 2016; available from https://github.com/laninsky/genetic_diversity_diffs). We included 12 sampled adult females and 12 adult males in this analysis, plus two females that could not be sampled but whose haplotypes were inferred from the haplotypes of their offspring (the adult female of Group 4, who supposedly had been replaced before the study period, and the adult female of Group 9, who was present during the study period but could not be sampled). Relatedness among females and among males was calculated using Wang’s estimator r and then compared using 1000 bootstrapping samples in Coancestry 1.0.1.9 91. In this analysis, as well as in the tests described below, we included 12 sampled adult females and 12 males. Individual heterozygosity was calculated using HL estimator (homozygosity by locus, see above) and compared between sexes using a paired t-test.
To evaluate spatial genetic structure, we conducted a spatial autocorrelation analysis following Smouse and Peakall 92 in PopGenReport 2.2.2 79, separately for adult females and males. The analysis calculated the correlation coefficient r between pairwise genetic distances, calculated using microsatellite genotypes with the method of Smouse and Peakall 92, and pairwise spatial distances, for each distance class. The coefficient r is bound between -1 and 1 and has a mean of zero when there is no correlation.  As a measure of spatial distances, we used distances between centroids of home ranges estimated using the 95% fixed kernel density method with ArcGIS Desktop 10.6 (ESRI). These distances varied from 215 to 3200 m.
To further evaluate spatial genetic structure in females, we conducted a test similar to spatial autocorrelation analysis using mtDNA haplotype distances, correlating the number of nucleotide differences between haplotypes with spatial distances. For this test, if a spatial genetic structure is present, a positive correlation between haplotype and spatial distances is expected. We used Mantel tests with 10,000 permutations in R package ecodist 93.
(6) Kinship within and between groups
Kinship was assigned conservatively based on r values and pedigree reconstruction in Colony. First-degree kinship (full-sibling and parent-offspring, without distinguishing between these two categories) was assigned to dyads with r > 0.487 (mean r for simulated parent-offspring dyads). Second-degree kinship (dyads sharing 25% of the genome, such as half-sisters or uncle-nephew dyads) was assigned for dyads with r > 0.247 (mean r for simulated half-offspring dyads). Dyads with lower r values were categorized as unrelated.