Computational protocol: Genetic Diversity and Population History of a Critically Endangered Primate, the Northern Muriqui (Brachyteles hypoxanthus)

Similar protocols

Protocol publication

[…] The sequences were aligned using CLUSTAL W and visually inspected in BIOEDIT . A BLAST search confirmed the absence of contamination with exogenous DNA. The B. hypoxanthus sequence from GenBank (AF213966) was used as a reference for alignment and haplotype determination. Haplotype determination and population comparisons were performed in ARLEQUIN 3.01 . Genetic variability was estimated by means of haplotype diversity (h) and nucleotide diversity (π) . The most appropriate model of nucleotide substitution was selected in jMODELTEST with the Akaike Information Criterion (HKY+I, base frequencies A = 0.3721, C = 0.2469, G = 0.1286, T = 0.2524, proportion of invariable sites of 0.9017, and transition/transversion ratio of 33). Whenever a model was required in subsequent analyses, either the best-fit model or the best approximation of it was used.Independent of current approximations of population sizes, we examined whether the number of haplotypes observed in large (2,000 ha or larger) versus small habitat patches was significantly different. We first inferred accumulation curves of haplotype richness for two sets of populations that were represented by at least 10 individuals. The first group was comprised of the populations found in areas larger than 2,000 ha (PERD, PESB, and SMJ). The SMJ population is comprised of small forest fragments interconnected with the total size of the area exceeding 2,000 ha (). The second group included the two populations found in areas of less than 1,000 ha (RPPN-FMA and RPPN-MS). To construct accumulation curves we generated 20 random samples of 10 individuals each for both of these two groups. Using these samples as input, we ran a rarefaction analysis in ESTIMATES 8.0 with seven cumulative sampling events of ten samples each and 1,000 randomization steps. The final seven samples, each consisting of ten individuals, summed to a simulated sampling effort of 70 individuals, which is approximately the largest combined N for each of these groups (first group N = 72, second group N = 75). The mean and standard deviations of the sampling intervals were plotted against the number of haplotypes with increasing sample size. A paired t-test was subsequently used to compare haplotype richness between these groups following the method employed by Peters et al. . [...] Analysis of Molecular Variance (AMOVA) was implemented in ARLEQUIN 3.01 using the model of Tamura and Nei to examine the distribution of genetic variance. AMOVA generates variance components and fixation indices that account for variation among populations (ΦCT), among social groups within populations (ΦSC), and within social groups among populations (ΦST). Statistical significance was assessed with 16,000 permutation steps. Populations represented by a small number of samples (PNC, FE and PEI) were excluded from this analysis. Because an exploratory phylogenetic analysis (not shown) indicated that there were no population sub-aggregations, the AMOVA included no division of groups of populations (the least inclusive hierarchy was the social group, and the most inclusive was the population).Because the RPPN-FMA and SMJ populations were sampled across different social groups, we tested whether both their overall and within-group haplotype compositions were randomly distributed. First, we compared the haplotype frequencies in the overall sample with a chi-square test. We assumed an extrinsic expectation in which all haplotypes should be found in equal frequencies proportionate to sample size (i.e. expected haplotype frequency  =  N haplotypes/N samples). To control for small sample size we first generated 10,000 random replicates from the sample frequencies (N = 64 for RPPN-FMA and N = 32 for SMJ) and then recalculated the chi-square for each replicate. The proportion of replicates with chi-square values greater than or equal to the observed value is the P-value of the test. We repeated this analysis for each social group assuming the expected frequencies of each haplotype to be equal to the observed population frequencies. In addition, to evaluate the contribution of group structure to the distribution of genetic variation we ran 16 AMOVAs, as described above, between these two populations. We removed one social group from each population at a time to detect groups that maximized the ΦSC. We used the serial goodness-of-fit plus method to correct for multiple testing ( evaluate whether social groups within the SMJ and RPPN-FMA populations are homogenous as a result of female dispersal, we used pairwise ΦST as a proxy for comparisons of between-group gene flow (see below). In addition to the four social groups previously defined for the RPPN-FMA population , we distinguished four social groups for the SMJ population by lumping individuals of different areas based on the groups’ proximity and the connectivity among forest patches (). The four SMJ groups defined were SSB, RP1/RP2, CO1/CO2, and RCT.Genetic differentiation due to isolation by distance was tested using the Mantel test , , which measures the relationship between genetic (pairwise ΦST) and geographic distances (Km) as implemented in IBDWS 3.15 . Statistical significance was calculated at 95% confidence interval (CI) with 10,000 randomization steps. We also ran the reduced major axis (RMA) regression analysis in IBDWS to graph the relationship between geographic and genetic distances. We inferred a haplotype network using the median-joining algorithm from NETWORK ( in order to depict haplotype relationships and explore the association of haplotype groups with geography. [...] We adopted a four-step approach to infer changes in the population size of northern muriquis. First, we calculated the neutrality tests FS , DF and F , and R2 , in DNASP 5.0 . These tests are suitable for detecting different types of demographic events (i.e. sudden expansion, sudden contraction, and bottleneck) in our data set . Very low and significant results of neutrality tests in mtDNA can be interpreted as evidence for population expansion. Significance was calculated after 50,000 coalescent simulations. Second, we calculated the exponential growth rate (g) in a coalescent framework using the Bayesian module in LAMARC 2.1.3 . One long Markov Chain Monte Carlo (MCMC) run was performed with 2 million generations, sampling every 20th step and a heating scheme of 4 simultaneous chains with different temperatures (1, 1.2, 1.5, and 3). Third, we used BEAST 1.4.8 to estimate the Bayesian skyline plot (BSP) to depict changes in effective population size over time. A normal distribution from the substitution rate estimated for Ateles (see Substitution rate) was used in BEAST. We ran the program for 100 million generations, with a sampling interval of 5,000 generations. Finally, we tested nine alternative scenarios for the diversification of four muriqui populations using the Approximate Bayesian Computation (ABC) approach implemented in DIYABC . We ran a total of 9 million simulations (1,000,000 per scenario) and used the logistic regression method to calculate the posterior probabilities of each scenario (see and for details). To minimize computation time and circumvent possible sampling bias that would affect the estimates, all coalescent analyses were carried out with a subset of 106 sequences. This smaller sample size was obtained after drawing 25 random samples from each RPPN-FMA and SMJ. Given the shallow haplotype network and lack of strong monophyletic subgroups (see ), all samples were treated as a single population in the first three treatments.The long-term female effective population size (Nef) was estimated using the formula Nef = θ/2μ, where μ is the substitution rate per site per generation and θ is the mutation-scaled effective population size. θ was estimated with the coalescent-based methods implemented both in LAMARC and MIGRATE-N 3.2.1 , , which utilized the maximum likelihood and Bayesian modules, respectively. The LAMARC analysis consisted of three independent MCMC runs of 10 initial chains of 30,000 steps each and 2 final chains of 2 million steps each sampling every 20th step. In MIGRATE, we ran 50 replicates of one long chain containing 10,000 samples each, sampling every 500th step. The prior distribution of θ was set as exponential with window (min. 0.0, mean 0.01, max. 0.1, delta 0.01). This analysis had a static heating strategy to explore the parameter space more effectively (swapping interval  = 1; four chains with start temperatures  = 1, 1.5, 3.0, and 1,000,000). The initial 10% genealogies were discarded as burn-in in all coalescent runs. Convergence for the coalescent analyses was assessed by inspecting the effective sample sizes and the stationarity of posterior distributions in TRACER 1.4.1 or the R statistical package . […]

Pipeline specifications