An Effort to Use Human-Based Exome Capture Methods to Analyze Chimpanzee and Macaque Exomes

Non-human primates have emerged as an important resource for the study of human disease and evolution. The characterization of genomic variation between and within non-human primate species could advance the development of genetically defined non-human primate disease models. However, non-human primate specific reagents that would expedite such research, such as exon-capture tools, are lacking. We evaluated the efficiency of using a human exome capture design for the selective enrichment of exonic regions of non-human primates. We compared the exon sequence recovery in nine chimpanzees, two crab-eating macaques and eight Japanese macaques. Over 91% of the target regions were captured in the non-human primate samples, although the specificity of the capture decreased as evolutionary divergence from humans increased. Both intra-specific and inter-specific DNA variants were identified; Sanger-based resequencing validated 85.4% of 41 randomly selected SNPs. Among the short indels identified, a majority (54.6%–77.3%) of the variants resulted in a change of 3 base pairs, consistent with expectations for a selection against frame shift mutations. Taken together, these findings indicate that use of a human design exon-capture array can provide efficient enrichment of non-human primate gene regions. Accordingly, use of the human exon-capture methods provides an attractive, cost-effective approach for the comparative analysis of non-human primate genomes, including gene-based DNA variant discovery.


Introduction
Non-human primates are increasingly studied as highly relevant animal models for human biomedical diseases and disorders. Members of the Macaca genus are among the most commonly studied non-human primates, due to their close evolutionary relationship to humans, analogous disease susceptibilities, and wide-spread commercial availability. The rhesus macaque (Macaca mulatta), estimated to have shared a common ancestor with humans approximately 25 million years ago (MYA) [1], is one of the most widely studied macaques. Genetic studies have shown the rhesus macaque to have common genetic risk factors with humans for age-related macular degeneration [2] behavioral disorders [3,4]_ENREF_4 and reproductive disorders such as amennorhea [5]. A close relative of the rhesus macaque, the Japanese macaque (M. fuscata) has served as a model for multiple sclerosis [6] and ischemia [7,8]_ENREF_8. The crab-eating or cynomolgus macaque (M. fascicularis) is widely used in studies of amyotrophic lateral sclerosis [9], and depression [10], among other disorders.
The chimpanzee (Pan troglodytes), is more closely related to humans than the macaques, sharing a common ancestor approximately 5-7 MYA [1]. The more recent divergence between humans and the chimpanzee has been of particular importance to the study of human evolution and speciation [11,12]_ENREF_12. In the field of comparative genomics, the chimpanzee genome provides a critical insight into studies of positive selection in primate genomes [13]. The chimpanzee has also served as a important model for neuroscience research, including studies of cognition [14], neurobiology [15], and behavior [16].
With the recent advance in genomic technologies, interest in comparative analysis of non-human primates, particularly as they relate to biomedical and evolutionary studies, has been rapidly expanding [17,18]. However, such studies are limited by the financial costs, computational requirements and effort required to generate genome-wide variant data on a large scale. Although improvements in next-generation sequencing (NGS) technology have already sharply reduced the cost of sequencing, the nonhuman primate still significantly lags behind in the comprehensive characterization of genome variation.
Exome sequencing has proven to be a powerful and efficient approach in human genetics studies [19], as it allows an unbiased investigation of almost all protein-coding regions in a large sample of individuals, at a fraction of the cost of whole genome sequencing. The method has been successfully applied to causative gene identification of several rare monogenic diseases such as Miller syndrome [20] spinocerebellar ataxias [21] and retinitis pigementosa [22]. A study of 50 Tibetan exomes uncovered a number of high-altitude adaptation related genes [23]. If the human exome-capture tools can be applied to the closely related non-human primate species, it could provide an opportunity to efficiently advance the pace of discovery of non-human primate sequence variants.
The human and chimpanzee genomes are about 99% identical, while macaques and human genomes are an estimated 93% conserved [17,18]. Given the high level of sequence conservation for coding regions among primates, we considered whether it would be feasible to efficiently enrich the exonic sequences of primate species using human-based exon capture designs. Applying exon-capture technology to non-human primate research would not only minimize cost, but it would also reduce the computational effort required for deep sequence analysis. Importantly, exome-sequencing approaches would expedite the discovery sequence variants of greatest interest to many investigators, those located in gene coding regions.
Similar efforts have been used to successfully enrich and sequence target regions of the Neanderthal genome [24]. More than a megabase of captured sequence was recovered from Neanderthal DNA, despite DNA degradation and the presence of significant microbial DNA contamination. This achievement provides support for the use of human exon-capture reagents for the study of more distantly related human ancestors.
Here we report an effort to use human based exome capture to analyze chimpanzee and macaque exomes. Nineteen non-human primates, involving 3 species, were evaluated. We report the utility of the human exon array tool for exon enrichment, DNA variant discovery, and for comparative genomic analysis.

Capture and Sequencing
We sequenced the exomes of nine chimpanzees (CM), two crabeating macaques (CE) and eight Japanese macaques (JP). Exonic sequences were enriched with the Agilent SureSelect all exon capture array (Human All Exon V1 for Human, CM and CE and Human All Exon V2 for JP)(Santa Clara, CA), targeting ,38 Mb (,46 Mb for JP) of DNA in nearly ,18,000 human consensus coding DNA sequences (CCDS). Sequencing was performed on an Illumina Hiseq2000 sequencer (San Diego, CA). The two human individuals (HM) were sequenced using the same workflow. The human genome [hg19, UCSC] was used as the reference for alignment of all sequenced individuals, since it enabled consistent comparisons for all individuals using the same coordinate system. All reads were aligned using SOAPaligner [25] with a gap-free model and for SNP calling and coverage calculations. BWA [26] gap tolerant alignments were also used to detect short indels. For the non-human primates, we also mapped reads to their own or nearest reference genome (chimpanzees to panTro2, crab-eating macaques and Japanese macaque to rheMac2) to evaluate sequencing quality (Table S2).

Evaluation and Comparison of Capture Performance
In order to perform unbiased evaluation, we compared the capture performance and sequencing quality in different species using equivalent metrics. Low quality reads and sequencing adaptor reads were filtered. Reads that mapped to the same chromosomal location and also had the same orientation were identified as duplicated reads. The one with highest mean quality score was retained and used in the analysis. The exomes of all 21 individuals were sequenced with a mean depth $28 fold (Table 1) on the array design target region (TR). Coverage of the TR ranged from 91.06% to 97.73%, with 67.59% to 81.33% of the sites in the TR covered by more than 10 reads. A gene-by-gene coding region coverage statistic was implemented for both theoretical TR and by actual sequence data. The theoretical analysis found that 83.68% (15559/18594) of CCDS genes' coding regions should have $90% covered by TR ( Table 2). The recovered sequence data showed that 80.40%-80.65% of the genes were covered by $90% of sequence reads for humans and chimpanzees, and by 77.76% for macaques ( Figure 1). Since coverage of chimpanzee CCDS was nearly equivalent to that of humans, we primarily focused on the macaques for our additional analysis. When mapping to the closest reference genome, the mapping rate was $84% in all 21 individuals, indicating a high quality sequencing data.
To further consider the exon sequence recovery, we used pairwise blastz alignment result of hg19/rheMac2 (download from UCSC genome browser) to evaluate the orthologous level between human and rhesus macaque in the TR. Evaluating a fragment from the target region of the human reference genome alignment to the Indian macaque reference genome, we calculated an orthologous score (OS) for every position and separated them into three levels: OS = 0, no alignment, unable to be captured theoretically; OS = 1, only one hit, high level of 1:1 orthologous coverage; OS = 2, multiple hits at different locations in the macaque reference, suggesting misalignment and potential for calling false polymorphisms (Table 3). We defined the region where OS = 1 as target orthologous region (TOR) for macaques and only used data in those regions in the following analysis. The region contains 93.5% of the initial target region. In TOR, 12,817 genes theoretically could be $90% covered (Table 2) and the actual sequence data showed that 88.86% (11389/12817) of them were $90% covered (Table S1). Thus based upon the TOR analysis, the coverage was close to the expected.
In assessing capture specificity, we defined the ratio of the reads aligned to a target region to the reads mapped to the species' closest reference genome. We found that capture specificity ranged from 40% to 74% (Table 1) and decreased as a function of sequence divergence. The percentage of clean reads that mapped to the TR for macaque were much less than that of human (,35% vs. ,68%).
In order to further investigate the influence of specific genomic features on capture efficiency, we evaluated how targets that failed to be captured differed from targets that were successfully captured (more than 50% of bases covered at least by one read) (Table 4A and 4B). We found that failed targets generally had more nucleotide differences from human (6.95% vs. 3.29%), higher percentage of indels (1.93% vs. 0.64%) and a higher GC content (57% vs. 47%) in macaques. More detailed inspection of GC content revealed that the total GC content influence capture rate over a broad range ( Figure 2 A and 2B). Target regions with moderate GC content (30%-50%) yielded higher coverage rates than regions with either high GC or low GC content ( Figure 3A). This finding is likely the result of the methods employed for the exon-capture, which included use of annealing temperatures that were optimize for the binding of target sequences with a moderate GC content.
Finally, we calculated the mismatch rate with the human and rhesus macaque reference genomes and correlated that with the average sequencing depth for all 158,852 autosomal targets. We found that sequencing depth decreased as sequence divergence increased ( Figure 3B).
The remaining unrecovered TOR in macaques totaled 1.74 Mb, accounted for by three potential factors: 1) there were ,0.22 Mb sites absent from both HM and CE in this capture array design; 2) the rhesus macaque genome was used to identify the TOR, and sequence divergence of the CE and JP genomes could have limited exon capture potential; 3) the amount of sequence, in combination with the reduced efficiency of capture, may have resulted in a biased representation of exon coverage. We conclude that we can predict the capture performance using available pairwise alignment information for those species that have their own/related reference genome. We observed TOR regions for Orangutan [ponAbe2] and Marmoset [calJac3] had a percentage of 92.26% and 88.80% separately in TR, using the same analysis with the IR genome. This demonstrates a potential utility of human based chip among these species. On the other hand for species without reference genomes, but which are very closely related to humans, exon capture should be sufficiently efficient to make sequencing the majority of the exome, as well as genotyping each individual, feasible.

Variation Discovery and Analysis
The use of exome sequencing for mutation discovery in human studies is critically dependent on the accurate identification of DNA polymorphisms and genotypes. We were interested in whether human exon capture technology would also enable the accurate detection of variants in non-human primates. Using the reference genomes of humans and non-human primates, it was possible to identify both variations between species and within species.
We used SOAPsnp [27] for calling genotypes. In order to directly compare different individuals, we only called genotypes in TORs where every individual had a read depth of $10X and had high quality reads. The qualifying regions totaled 17.07 Mb, or 45.36% of the original TR (Table 2). In total, 79,704 ,311,477 single base variants were detected within non-human primates (compared with human reference) and 15,013,27,391 of them were candidate intra-species polymorphisms (Table 5A).
We considered the possibility of cross DNA contamination contributing to the SNPs called. Though we used care in the DNA extraction, capture, library construction and sequencing procedures, we nonetheless evaluated whether the derived exome Sequencing data showed that 80.40%-80.45% of genes were covered $90% for both humans and chimpanzees, and 77.76% for macaques. Coverage of macaque CCDSs was higher than target orthologous region. The coverage of human and chimpanzee genes was close when concerns target region, and higher for macaque compared to target orthologous region. doi:10.1371/journal.pone.0040637.g001 sequence contained evidence of human DNA contamination. We compared all exome sequencing data to the closest reference sequence. We found that the exome sequences and closest reference sequences exceeded a 99% match in all 19 non-human primate individuals (Table 5A). These results indicate that the sequence data used for genotype calling aligned most closely with the reference genome and did not suggest the presence of contaminating human DNA. The abundance of polymorphisms discovered in coding and flanking regions enabled us to explore the potential evolutionary history of these species on a population scale. We found the macaques had the highest level of heterozygosity, followed by the chimpanzee. These data suggest that the macaques have a higher effective population size than chimpanzees, which in turn have a higher effective size than humans (Table 5B). In all species the ratio of nonsynonymous to synonymous SNPs was larger than the ratio of nonsynonymous to synonymous differences (Neutrality index, Table 5B) This is consistent with previous studies of effective population sizes in these species, [13,17] suggesting that nonsynonymous deleterious alleles are the subject of selection pressure within each species.
Since the nineteen non-human primate exomes included three different species, we used the polymorphism data to evaluate the evolutionary relationships among individuals. We constructed a phylogenetic tree (Figure 4) using the neighbor-joining (NJ) method. The phylogenetic tree is consistent with the ancestry analysis shown in previous genome studies [17,18], as well as with results based upon microsatellite markers [1]. Principal component analysis (PCA) also showed that the individuals in each species were genetically most similar ( Figure 5).
We also identified short insertions and deletions (indels) in the exons, which are likely to be functionally important and may contribute to species divergence (Table S3). Here we mapped the reads from each individual to the human reference sequence using BWA and Samtools [28] to identify indels (see methods for details). In total, 140,1,883 (range for all individuals) human-primate coding indels were discovered in each individual (Table 6), and 54.61%,77.32% of them were 3 base pairs in length. The distribution of indel lengths in the coding regions is consistent with a previous study [29] across all of the 21 exomes (Figure 6), reflecting the evolutionary pressure to preserve intact reading frames.
We carried out Sanger sequencing for validation on a randomly selected 41 SNPs and 18 indels in two exome sequenced CE individuals. The validation rate was 85.4% for SNPs and 84.2% for indels (Table S4), indicating a high accuracy rate.

Discussion
Here we evaluated the efficiency of exome capture and sequencing in non-human primates using a human capture array. We determined that although capture specificity decreased as sequence divergence increased, it is still a viable option for crabeating macaque, Japanese macaque and Chimpanzee exome enrichment. We also identified inter-species and intra-species variation by comparing recovered sequences to the rhesus macaque reference genome, and then validated the call accuracy using Sanger resequencing.
One limitation of this method is the inability to detecting large structural variations, such as large indels, inversions and copy number variations (CNV). In addition, the particular human exon capture design may miss some important 59 and 39 noncoding as well as small RNAs. Further, a subset of the genomic target regions were not recovered.
Despite the limitations, use of human exome-enrichment has some distinct advantages. We were able to efficiently capture gene coding region for both species with reference genomes, as well as those without a species-specific reference genome. Applications of this method could enable the efficient, large scale exon-sequencing of non-human primates, potentially identifying DNA variants of relevance to human disease, and guiding animal model development and translational research. Finally, exome-capture approaches could expedite the comparative genomic analysis of non-human primate species, for example between gorilla, orangutan, gibbons, baboon, and even New World monkeys, providing insight to the evolution of the human genome.
The nineteen non-human primate exomes presented here also highlight the degree of variation existing between these widely used non-human primate animal models. The abundant genetic diversity evident in individual primates from distinct geographic populations is of direct interest to primatology, medical research, population genetics and phylogeographic studies.

Ethics Statement
Blood samples of 2 Han Chinese individuals were collected from Anhui Medical University. We have obtained ethics approval for our study from the Ethical Committee of Anhui Medical University. We have also obtained informed written consent from all participants involved in this study. The authors of this paper did not collect the blood samples themselves; all non-human primate blood samples used in this study were initially collected for routine health check purposes.
Samples from nine wild born, unrelated chimpanzees (Pan troglodytes schweinfurthii) were collected at the Ngamba Island Chimpanzee Sanctuary in Uganda. The samples were taken during routine health checks by the sanctuary veterinarians and were exported under CITES export permit Sn.UG 002249. The Ugandan Wildlife Authorities and Uganda National Council for Science and Technology approved the research that resulted in the CITES export permission. Figure 2. GC distribution of captured and un-captured targets by CE1. Sequenced CE1 exome was compared with its closest reference genome (rheMac2) to examine the influence of GC content on coverage. (A) Captured targets were defined as targets with more than 50% of bases covered at least by one read. The captured targets were mainly of moderate GC content. (B) Uncaptured targets were less than a half of the bases coverd by at least one read. The uncaptured targets were dominated by high GC content, except for a small portion with low GC content. doi:10.1371/journal.pone.0040637.g002 Figure 3. Average sequencing depth of captured targets versus human reference and rhesus macaque reference mismatch. Mismatch rate of each target region was calculated by comparing human reference and rhesus macaque references respectively. (A) The X-axis identifies the mismatch rate region, and the Y-axis identifies the mean depth of target region, as calculated from CE2. As the mismatch rate increases, the mean depth decreases. (B) Average sequencing depth of captured targets versus each target GC content. Targets with extreme GC content were poorly captured; the mean depth was lower compared with moderate GC content. doi:10.1371/journal.pone.0040637.g003 Table 5.  China) under the temperature of, between 16uC and 29uC and 50,80% relative air humidity, with four air changes per hour and a 12 h light:12 h dark cycle. Animals were fed apples (100 g/day each) and monkey Chow (Feed Research Institute, Guangzhou, Guangdong) twice daily and allowed free access to water. The blood samples were collected during routine veterinary health checks. Guangdong Landau Biotechnology Co. is accredited by the Association for Assessment and Accreditation of Laboratory Animal Care International (AAALAC). The macaques were originally exported from Kampuchea under CITES export permit Sn. IC 0381. All macaque experiments were subject to approval and surveillance by the Institutional Animal Care and Use Committee of Guangdong Landau Biotechnology Co., Ltd. Genomic DNA was extracted from blood samples using the Nucleon kit (TaKaRa, Japan). Blood derived DNA was used to minimize somatic and cell-line derived false positives.
Eight Japanese macaques (Macaca fuscata) were born and raised at the Oregon National Primate Research Center (ONPRC); all animal procedures were approved by the ONPRC Institutional Animal Care and Use. Committee and conformed to the NIH guidelines on the ethical care and use of animals in research. The selected animals were members of captive breeding group that lived within in a 2 acre, outdoor corral, and consumed Primate Diet no. 5000, (Lab Diet, Richmond, IN). Blood samples were collected during routine veterinary care, by venipuncture and collection into an EDTA vacutainer. Genomic DNA was isolated using the Wizard DNA Extraction Kit (Promega, Inc.).

Exon Capture and Sequencing
The qualified genomic DNA samples were randomly fragmented by Covaris and the sizes of the library fragments were distributed between 150 to 200 bp. Adapters were then ligated to both ends of the resulting fragments. The adapter-ligated templates were purified by the Agencourt AMPure SPRI beads and fragments with insert size of about 250 bp were excised. Extracted DNA was amplified by ligation-mediated PCR (LM-PCR), purified, and hybridized to the SureSelect Biotinylated RNA Library (BAITS) for enrichment. Hybridized fragments were bound to the strepavidin beads whereas non-hybridized fragments were washed out after 24 h. Captured LM-PCR products were subjected to Agilent 2100 Bioanalyzer to estimate the magnitude of enrichment.
Each captured library was loaded on the Illumina Hiseq2000 platform for high-throughput sequencing to the desired average sequencing depth (The exomes of all 21 individuals were sequenced with a mean depth $28 fold on the array design target region). Raw image files were processed by Illumina basecalling Software 1.7 for base-calling with default parameters and the sequences of each individual were generated as 90 bp pair-end reads. Read Mapping and Genotype Calling SOAPaligner (soap 2.21) was used to align reads to the human reference genome (hg19) allowing a maximum of 3 mismatches per 90 bp read. Full SOAP parameters were: -a -b -D -o -2 -t -v 3 -l 35 -s 40 -m 0 -x 500 -p 4 -r 1 (-v 5 for macaques). Reads that aligned to the designed target region (TR) were collected for genotype calling and subsequent analysis.
Based on SOAP alignment results, the software SOAPsnp was used to call genotypes. The following parameters were set: -r 0.0005e 0.001 -t -u -2 -i -d -o -M -L 90 -s -T(http://soap.genomics.org.cn/ for details).
Raw data has NCBI Short Read Archive accession no. SRA038809.

SNP Identification
Single base differences between human and non-human primates references (
Alignment of fragments from the TR between human and nonhuman primates were used to calculate OS. Alignment hits count = 0, OS = 0; alignment hits count = 1, OS = 1; alignment hits count$2, OS = 2. Regions in TR where OS = 1 were defined as TOR.

Phylogenetic Tree
We use high quality genotypes on chromosome 10 (TOR, depth$10 and consensus quality$20 in every individual) generated by SOAPsnp to build a phylogenetic tree with the program

Identification of Insertions and Deletions
Gap tolerant alignments to human reference (hg19) were used to call indels with program BWA using parameters: aln -o 1 -e 63 -i 15 -L -l 31 -k 2 -t 4. Insertions and deletions (indels) were identified using samtools, with the command lines as following: samtools mpileup -ugf ref.fa -b bam.list | bcftools view -bvcg -. var.raw.bcf bcftools view var.raw.bcf | vcfutils.pl varFilter -D10000. var.flt.vcf. Filter criteria were $3 reads support and number of indel supported reads $30% of all reads mapped to the genomic position. Indels were called as heterozygous if the indel supported reads were 30-70% of all reads at that position, and homozygous if they were greater than 70%.

Validation of SNP and Indel Variants
Forty one SNPs and 18 indels were selected randomly for validation of SNPs and indels. The selected SNP and indels were genotyped by PCR and Sanger sequencing. Primers for each selected SNP and indel were designed based on the IR genome; the detail sequences of each primer pairs were supplied in Table  S4. The polymerase chain reaction (PCR) was performed in a final volume of 50 ul with 30 cycles at 94uC for 30 s, annealing temperature (Table S4) for 30 s, and 72uC for 30 s. The PCR product was purified and sequenced by BGI (BGI-Shenzhen, Shenzhen 518083, China).