|Home | About | Journals | Submit | Contact Us | Français|
Conceived and designed the experiments: XJ HD Xiaoning Wang Jian Wang Jun Wang HZ YW. Performed the experiments: BF FS L. Sun MZ Jufang Wang FL XL FJ Xinzhong Wang BL Y. Zhang JH Jing Wang XZ MHS. Analyzed the data: XJ MH BF YM LO JR JS L. Song Y. Zhu XF GZ. Wrote the paper: XJ BF TM MHS HD. Prepared samples: BF FS L. Sun MZ Jufang Wang FL XL FJ Xinzhong Wang BL Y. Zhang JH Jing Wang XZ MH. Assisted in writing manuscript: XJ MH BF YM LO JR TM FS L. Sun JS MZ L. Song Jufang Wang FL Y. Zhu CH HS XL ZG FJ Xinzhong Wang BL Yu Zhang JH Jing Wang HZ YW XF GZ Jian Wang XZ MHS HD Jun Wang Xiaoning Wang.
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.
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) , 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  behavioral disorders , _ENREF_4 and reproductive disorders such as amennorhea . A close relative of the rhesus macaque, the Japanese macaque (M. fuscata) has served as a model for multiple sclerosis  and ischemia , _ENREF_8. The crab-eating or cynomolgus macaque (M. fascicularis) is widely used in studies of amyotrophic lateral sclerosis , and depression , among other disorders.
The chimpanzee (Pan troglodytes), is more closely related to humans than the macaques, sharing a common ancestor approximately 5–7 MYA . The more recent divergence between humans and the chimpanzee has been of particular importance to the study of human evolution and speciation , _ENREF_12. In the field of comparative genomics, the chimpanzee genome provides a critical insight into studies of positive selection in primate genomes . The chimpanzee has also served as a important model for neuroscience research, including studies of cognition , neurobiology , and behavior .
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 , . 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 non-human 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 , 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  spinocerebellar ataxias  and retinitis pigementosa . A study of 50 Tibetan exomes uncovered a number of high-altitude adaptation related genes . 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 , . 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 . 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.
We sequenced the exomes of nine chimpanzees (CM), two crab-eating 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  with a gap-free model and for SNP calling and coverage calculations. BWA  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).
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 11 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.
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  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 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, ,  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 , , as well as with results based upon microsatellite markers . 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  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  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.
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 crab-eating 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 5′ and 3′ 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.
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.
Two adult cynomolgus macaque (Macaca fascicularis) were individually housed in cages with the size of 70 cm×70 cm×80 cm by the South-China Primate Research & Development Center (Guangdong Landau Biotechnology Co., China) under the temperature of, between 16°C and 29°C 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.).
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.
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.0005 -e 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.
Single base differences between human and non-human primates references (ref-ref variants) were identified from pairwise blastz alignments result. *.axt files for Human/Chimp (panTro2) and Human/Rhesus (rheMac2) were downloaded from UCSC Genome Browser (http://hgdownload.cse.ucsc.edu/downloads.html#human).
Raw point variants from exome sequencing data of each individual were extract from SOAPsnp generated genotype with filter criteria Q20 and depth ≥10. Homozygous variants consistent with ref-ref variants were defined as interspecies variants. Remaining raw point variants were considered to be potential intraspecies polymorphisms.
Orthologous score (OS) for every site in TR was calculated with pairwise blastz alignments result. The *.axt files for Human/Chimp, Chimp/Human (panTro2) and Human/Rhesus, Rhesus/Human (rheMac2) were downloaded from UCSC Genome Browser (http://hgdownload.cse.ucsc.edu/downloads.html#human).
Alignment of fragments from the TR between human and non-human 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.
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 TreeBeST (treebest-1.9.2, http://treesoft.sourceforge.net/treebest.shtml; L.Heng, A.J. Vilella, E.Birney, and R.Durbin, in prep) with the model of neighbour-joining tree, SDI, rooting (threshold: nj -b 500).
We use high quality genotypes on chromosome 10 (TOR, depth≥10 and individuals whose quality score under 20 were discarded) generated by SOAPsnp to perform principal component analysis (PCA) using EIGENSOFT3.0 (parameter: -i all.evec -c 12 -p HM:CM:CE:JP -x).
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%.
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 94°C for 30 s, annealing temperature (Table S4) for 30 s, and 72°C for 30 s. The PCR product was purified and sequenced by BGI (BGI-Shenzhen, Shenzhen 518083, China).
Coverage of targeted orthologous genes. Coverage of 18,594 gene coding regions of human CCDS by theoretical target region, target orthologous region and sequencing data from 3 species.
Data production. Summary of captured target sequence coverage for each non-human primate exome and two human exomes.
Indels. Indels of 21 individuals identified through aligning the exomes to the human reference genome.
Sanger sequencing validation result of SNP and Indels.
We thank Prof. Rasmus Nielsen from UC-Berkley for his suggestions and comments. We are also indebted to many additional faculty and staff of BGI-Shenzhen who contributed to this teamwork. For generously providing chimpanzee samples, we thank Ngamba Island Chimpanzee Sanctuary.