The common marmoset genome provides insight into primate biology and evolution.

We report the whole-genome sequence of the common marmoset (Callithrix jacchus). The 2.26-Gb genome of a female marmoset was assembled using Sanger read data (6×) and a whole-genome shotgun strategy. A first analysis has permitted comparison with the genomes of apes and Old World monkeys and the identification of specific features that might contribute to the unique biology of this diminutive primate, including genetic changes that may influence body size, frequent twinning and chimerism. We observed positive selection in growth hormone/insulin-like growth factor genes (growth pathways), respiratory complex I genes (metabolic pathways), and genes encoding immunobiological factors and proteases (reproductive and immunity pathways). In addition, both protein-coding and microRNA genes related to reproduction exhibited evidence of rapid sequence evolution. This genome sequence for a New World monkey enables increased power for comparative analyses among available primate genomes and facilitates biomedical research application.

Apparently unique among mammals, marmosets routinely produce dizygotic twins that exchange hematopoietic stem cells in utero, a process that leads to lifelong chimerism[1,2]. As a result of this placental exchange, the blood of adult marmosets normally contains a substantial proportion of leukocytes that are not derived from the inherited germline of the sampled individual but rather were acquired in utero from its co-twin. In addition, marmosets (subfamily Callitrichinae) and other callitrichines are small in body size as a result of natural selection for miniaturization. This reduced body size may be related to their gestation of multiples and to their social system, also unique among primates[3-5]. These animals display a cooperative breeding system in which generally only one pair of adults in any social group are active breeders. Other adult group members participate in the care and feeding of infants, but do not reproduce. This alloparental care is rare among anthropoid primates, with the clear exception of humans. The evolutionary appearance of major new groups (e.g. superfamilies) of primates has generally been characterized by progressive increases in body size and lifespan, reductions in overall reproductive rate, and increases in maternal investment in the rearing of individual offspring. In contrast, marmosets and their callitrichine relatives have undergone secondary reduction in body size from a larger platyrrhine ancestor[6], and have evolved a reproductive and social system in which the dominant male and female monopolize breeding but benefit from alloparental care provided to their offspring by multiple group members.

Here we report an analysis of the complete genome sequence for this first New World monkey (NWM) to be sequenced (Supplementary Note). Our results include comparisons of this platyrrhine genome with the available catarrhine (human, other hominoid and Old World monkey (OWM)) genomes that reveal previously undetected aspects of catarrhine genome evolution, including positive selection in specific genes and significant conservation of previously unidentified segments of noncoding DNA. The marmoset genome displays a number of unique features, such as rapid changes in microRNAs (miRNAs) expressed in placenta, and non-synonymous changes in protein-coding genes involved in reproductive physiology, which may be related to the frequent twinning and/or chimerism.

The genetics of twinning

WFIKKN1, a multidomain protease inhibitor that binds growth factors and BMPs[7], has non-synonymous changes found exclusively in common marmosets and all other tested callitrichine species that twin. In the one callitrichine species that does not produce twins (Callimico), one change has reverted back to the ancestral (non-twinning) primate sequence. GDF9 and BMP15, genes associated with twinning in sheep and humans, also exhibit non-synonymous changes in callitrichines.

Genetics of body size

We detected positive selection in five growth hormone/insulin-like growth factor (GH-IGF) axis genes with potential roles in diminutive body size and eight genes in the nuclear-encoded subunits of the respiratory complex I that affect metabolic rates and body temperature, adaptations associated with challenges of small body size.

Novel microRNA loci

Marmosets exhibit a significant number of unanticipated differences in miRNAs and their targets, including 321 novel miRNA loci. Two large clusters of miRNAs expressed in placenta show substantial sequence divergence compared with other primates and are potentially involved in marmoset reproductive traits. We identified considerable evolutionary change in the protein coding genes targeted by the highly conserved let-7 family and significant co-evolution of the rapidly evolving Chr. 22 miRNA cluster and their targets.

Novel functional elements specific to primates

The marmoset genome provides unprecedented statistical power to identify sequence constraint among primates, facilitating discovery of novel genomic regions underlying primate phenotypic evolution. The 23,849 regions that exhibit significant sequence constraint among primates but not in non-primate mammals are overwhelmingly noncoding, disproportionately associated with genes involved in neurodevelopment and retroviral suppression, and frequently overlap transposable elements.

Novel detection of positive selection in human, ape and OWM lineages

For seven genes, we detected positive selection on the branch leading to Catarrhini. Five were previously unidentified, including genes involved in immunobiology and reproduction (Table 1).

Table 1

Gene Ontology (GO) Categories Enriched for Positively Selected Genes.

		Genes			Adjusted
GO(a)	Description	PSG(b)	Total(c)	Excess(d)	P-value (MWU)	P-value (Holm)	P-value (FET)
0005576	extracellular region	150	1954	1.3	3.24E-15	9.80E-12	3.86E-17
0005615	extracellular space	63	429	2.4	2.52E-08	7.61E-05	1.31E-08
0005747	mitochondrial respiratory chain complex I	8	14	9.4	1.81E-07	5.47E-04	2.72E-05
0006952	defense response	54	324	2.7	2.19E-06	6.59E-03	3.38E-09
0004872	receptor activity	103	866	2.0	3.42E-06	1.03E-02	1.05E-08
0007606	sensory perception of chemical stimulus	20	136	2.4	5.82E-06	1.75E-02	1.26E-03
0030246	carbohydrate binding	29	203	2.3	6.81E-06	2.05E-02	1.78E-04
0006954	inflammatory response	36	181	3.3	8.39E-06	2.52E-02	3.31E-08
0004984	olfactory receptor activity	16	107	2.5	9.88E-06	2.97E-02	3.21E-03
0009611	response to wounding	53	332	2.6	2.93E-05	8.79E-02	1.73E-08
0006955	immune response	41	295	2.3	3.18E-05	9.53E-02	1.57E-05

GO category number,

positively selected genes (PSG) identified with a threshold of P<0.05,

total genes in the GO category,

fold enrichment of PSGs over background. Enriched GO categories were identified by Mann-Whitney U-test (MWU), nominal P-value adjusted for multiple testing by Holm correction (Holm), and Fisher’s exact test (FET) using all genes with nominal P<0.05 in marmoset lineage likelihood ratio test. Note that the results of MWU may also be affected by relaxation of constraint, while FET considers only genes identified as under positive selection.

Results

Genome Assembly and Features

The 2.26 Gb genome of a female marmoset (#186/17066) assembled with Sanger read data (6x) and a whole genome shotgun (WGS) strategy (Supplementary Tables 1-4, Figure 1) represents ~90% of the marmoset genome. By all avilable measures the chromosomal sequences have high nucleotide and structural accuracy (Contig N50 29 kb, Scaffold N50 6.7 Mb, Supplementary Notes) and provide a suitable template for initial analysis.

Given the inherent genetic chimerism in this species, blood DNA contains sequences from the germline of the sampled individual and also her male co-twin. We took advantage of the sex-difference in the co-twins to estimate the proportion of reads originating from the co-twin (see Supplementary Tables 5-6, Figure 2, Note). These analyses indicate that 10% of the reads in the reference genome dataset are derived from the co-twin.

We estimated the amount and size of marmoset segmental duplications using two computational methods, WGAC[8] and WSSD[9]. The assembly based duplications add a total of 138 Mb of nonredundant sequences (4.7% of the whole genome), slightly less than observed in human or chimpanzee (~5%)[10-12] but more than orangutan (3.8%)[13] where specific collapses in the released assembly version might explain this anomaly. See Supplementary Figure 3, 4, Tables 7-10 and Note.

For segmental duplications >10 kb with >94% sequence identity (Supplementary Table 8) we compared the two independent methods to measure artifactual duplications and mistaken assembly collapses. Both methods identified a total of 18 Mb of duplications, while 26 Mb were possible artifactual duplications and 53 Mb were possible collapses. To validate the methods we tested 97 clones by FISH mapping to marmoset chromosomes (see Supplementary Table 9). Bothe methods successfully identify segmentally duplicated regions and unlike previous Table studies, WGAC seems better suited than WSSD to detect duplication in marmoset. The degree to which this is due to the chimeric nature of the individual sequenced is not clear, though the chimerism is certainly a contributing factor.

The overall repeat composition of the marmoset genome is similar to other sequenced primate genomes[10,12-14], containing ~1.1 million Alu elements, ~ 660,000 of which are full-length. However, in the recent past, Alu retrotransposition appears to be somewhat slower in marmosets than in humans and rhesus macaques (Supplementary Note).

Constrained sequence evolution indicates natural selection, and therefore implies conserved function. By extension, lineage-specific constraint indicates lineage-specific function[15,16]. Using the marmoset we detected 23,849 elements constrained in anthropoid primates but not in non-primate mammals[17] (Supplementary Note). These anthropoid-specific constrained (ASC) sequences potentially drove primate phenotypic evolution, and are abundant in non-coding regions (e.g. upstream of SNTG1) (Supplementary Figures 5a, 5b) though coding exons are also represented (e.g. PGBD3). Annotated transposable elements contributed 46% of ASC basepairs. We validated the enhancer activity of six (of eight tested) elements in human embryonic stem cells (Supplementary Figure 5c, Table 11), and showed that their mouse orthologs had little or no functional activity. This dataset highlights specific loci that acquired new functional roles in the primate lineage, and potentially reveals molecular mechanisms underlying unique primate traits.

Gene Content and Gene Families

The Ensembl gene set[18] (Supplementary Note) of 21,168 genes (44,973 transcripts) included 219 genes with marmoset protein support and 15,706 genes without marmoset protein evidence but with human protein evidence. The remaining 5,243 genes had transcripts supported by protein data from other sources (Supplementary Figures 6g, 6h).

A phylogenetic framework including four other primates, two rodents and three Laurasiatheria, showed 429 primate-specific gene families, among which few were present only in marmoset (Supplementary Figure 7, Tables 12-19, Note). More than half of these families (221/429) are indeed absent in marmoset, suggesting that they emerged after the catarrhine-platyrrhine divergence. In addition, many families are absent in rhesus macaque, and thus almost half are apparently unique to apes.

Our comparative analysis found surprising changes in the miRNA repertoir and the mRNA targets that they regulate. We identified 777 mature miRNAs (mapping to 1,165 hairpin precursor miRNAs) (Supplementary Table 20). Most were confirmed through expression studies (582, 75%) (Supplementary Note) and conserved in primates (~55-58%). Many (321 miRNAs mapping to 477 hairpins) are novel (not found in any other species to date). These include miRNAs exclusive to marmoset, miRNAs exclusive to Platyrrhini, and conserved miRNAs that are yet to be discovered in other species. The two largest marmoset miRNA clusters (Chr. 22, Chr. X) are expanded in number compared to human (112 vs. 49 human Chr. 22 hairpins and 40 vs. 15 human Chr. X hairpins) (Supplementary Table 22) and diverged in sequence. Less than 3% of Chr.22 and 8% of the Chr. X miRNAs are conserved across primate (Supplementary Table 22), and most exhibit at least one nucleotide modification in the 5′-seed region (83% of Chr. 22, 78% of Chr. X) compared to human (Supplemental Tables 20, 22, 23, 29). The rapidly evolving Chr. 22 and Chr. X clusters dominate the miRNA expression in marmoset placenta whereas the marmoset brain exhibits a more diverse miRNA expression pattern (Supplementary Figure 9, Tables 30-32). In contrast, some families (e.g. let-7) are completely conserved in all five primates (Supplementary Figure 8).

Changes in the miRNA seed region are expected to correspond to changes in the genes they regulate unless the miRNAs and their mRNA targets have co-evolved. Comparing the annotated genes containing predicted let-7 target sequences (Figure 1, Supplementary Note), we found 165 common to human and marmoset, 44 unique to marmoset and 64 unique to human. Despite caveats related to differences in assembly and annotation qualities, it is striking that less than half of the targets for this highly conserved family are shared between marmoset and human (Supplemental Table 34), a number similar to non-euarchontoglires (dog, horse and cow). A phylogenetic analysis of these changes shows that let-7 targets have evolved rapidly in primates compared to other species (Figure 2). The pattern of miRNA:mRNA target evolution differs among the three described miRNA families and even between the two rapidly evolving families (Supplementary Tables 33-37). In the Chr. X cluster, as expected, fewer than 50% of the target sequences are shared between marmosets and humans (Supplementary Table 35). In contrast, in the Chr. 22 cluster 84% of the targets are shared (Supplementary Table 36), implying significant co-evolution of miRNAs and their targets in the Chr. 22 cluster but not in the Chr. X cluster.

Figure 1

Predicted let-7 Regulated Genes (miRNA targets)

The number of protein coding genes with predicted targets for let-7 miRNA binding in the 3′ UTR are shown. Only single copy orthologs are counted and numbers are relative to the number found in human (100% on scale). The number of gene targets shared with human (in blue) falls as the evolutionary distance increases, as expected. However, the proportion of let-7 targets shared with human is comparable for marmoset, dog, horse and cow, while mouse and rat share fewer targets with human than other non-primate placental mammals.

Figure 2

Gains and Losses of let-7 Regulated Genes

The conserved let-7 miRNA targets variable numbers of genes. Let-7 targets gains (in green) and losses (in blue) mapped to the phylogenetic tree of the analyzed species, line thickness indicates the rate of gain or loss. Gains and losses that occurred twice on independent lineages were omitted. Gains exceed losses on each branch of the tree, and the total gained (196) is four times the losses (49). Primate lineage changes (gains plus losses) exceed non-primate lineage changes (except for the branch leading to rat after divergence from mouse).

Small marmosets are believed to have evolved from a larger ancestor, so we looked for positively selected genes that might explain this change. We identified 37 positively selected genes on the marmoset lineage and 7 on the branch to Catarrhini (FDR<0.01)(Supplemental Table 38). Five of these seven genes (SAMHD1, CLEC4A, ANKZF1, KRT8, CATSPERG) were previously unrecognized as positively selected[19]. An additional 91 positively selected genes could not be traced to a particular branch due to a lack of identifiable outgroup orthologs. Following trends observed in previous studies[19], gene ontology categories related to immunity, physiological defense response, and sensory perception are enriched (Table 1). In addition, the ATP synthesis and transport and NADH dehydrogenase activity categories show enrichment (Mann-Whitney U-test, P<0.05). The latter group contains eight positively selected nuclear genes encoding subunits of the respiratory Complex I. Resulting differences in Complex I regulatory and kinetic properties could affect metabolic rates and body temperature, challenges posed by small body size.

A prominent example of marmoset lineage positive selection can be found in IGF1R (P=0.0014), which is associated with short stature in humans[20,21]. This protein shows multiple mutations within crucial binding domains (Figure 3) which likely affect ligand-receptor binding affinity. Other growth hormone related positively selected genes possibly related to small stature include growth hormone secretagogue receptor (GHSR), insulin-like growth factor 2 (IGF2), insulin-like growth factor binding proteins 2 (IGFBP2), and 7 (IGFBP7), and epidermal growth factor (EGF) (marmoset lineage, P<0.05). Targeted exon sequencing of multiple species revealed several callitrichid-specific non-synonymous substitutions in genes that are strong candidates for influencing diminutive body size (GDF9, BMP15, BMP4). Analysis of these mutations by SIFT[22] and Polyphen[23] indicate that these likely affect the function of the corresponding proteins[24] (Supplementary Note, Table 38).

Figure 3

Residues under Positive Selection in IGF1R

The insulin-like growth factor 1 receptor (IGF1R) interacts with other genes in growth hormone pathways and has a role in both prenatal (left) and postnatal (right) growth. Genes in these pathways in marmoset that have residues under positive selection are tallied, the number changes that can be assigned to either the marmoset or callitrichine NWM lineages are also shown. In the middle, the first three domains of the IGF1R alpha chain are shown with positively selected residues in red (Bayes empirical Bayes analysis PP>0.95) and yellow (PP>0.5). Leucine-rich repeat domains L1 and L2 are shown in green with L1 on top, cystein-rich region CR is shown in blue. A multiple alignment of IGF1R proteins from several mammalians species (bottom) exhibits several marmoset changes in a short region corresponding to the part of structure enclosed in the black rectangle.

The genetic basis of twinning has received significant attention in humans and other animals[25-27]. Genetic differences drive variation in ovulation number among sheep strains[25,28]. There is also clear evidence for genetic influence on human twinning but specific genes have not been identified. We studied 63 candidate genes previously implicated in the control of either body size, number of ova produced in a single estrous cycle, or both. Of these, 41 genes with putative marmoset-specific non-synonymous variants were examined further (Supplemental Table 39, 40). Three genes with a role in ovulation (BMP4, FSTL4 and WFIKKN1) showed likely function-altering amino acid changes as scored by both SIFT[22] and Polyphen[23] (Supplementary Note and[24]). Potentially functional non-synonymous substitutions in genes coding for follicle stimulating hormone receptor (FSHR), BMP10, BMP15, GDF9 and GDF15 were also found. Notably, a single non-synonymous substitution in WFIKKN1 is common to all callitrichids we tested, with the exception of Callimico goeldi. That species shows a reversal of this change back to the sequence found in OWM and other non-twinning NWM. C. goeldi is the only callitrichid that does not regularly twin, and given its phylogenetic position, it is highly likely to have reverted back to singleton births from an ancestral state that exhibited twinning. This amino acid change in WFIKKN1 is therefore a strong candidate for a role in the origin of twinning in callitrichids.

Hematopoietic chimerism of marmosets was expected to correlate with significant changes in immune system function. We found positively selected genes related to immune response significantly enriched in marmoset (Table 1). NAIP and NLRC4 homologues, conserved in mammals, are absent in marmoset (Supplemental Table 38). These proteins form the NAIP inflammasome in macrophages, a cytoplasmic complex that triggers macrophage inflammatory death through activation of caspase-1[29,30] and could affect reproduction as human NAIP is expressed in the placenta.

Other positively selected genes potentially involved in circumventing unwanted chimerism-associated responses include CD48, a ligand for CD244 (2B4) found on the surface of hematopoietic cells which regulates natural killer cells[31] and interleukins IL5 and IL12B, involved in T cell development and in allergic responses[32]. Finally, in contrast to the extensive family of KIR genes integral to immune system function in humans and other catarrhine primates, the marmoset genome contains only two KIR genes, one of which is partial.

Most differences in protease gene families observed between marmosets and other primates occur in genes related to reproductive and immune systems (Supplementary Note). For example, ADAM6, with a role in fertility[33,34], was lost in marmoset, while ISP2, involved in embryo implantation[35], has been duplicated twice. KLK2/3, duplicated in the Catarrhine ancestor[36] and involved in reproductive physiology[33] is non-functional in marmoset. Chymase and tryptase protease changes and CMA1 and MAST duplications potentially affect the immune response[37,38] and mast cell biology, respectively. The duplicated CMA1 gene may be related to the murine-specific mast cell proteases (MCPs) absent in hominoids[39]. Changes in the C-terminal of MMP19 an IGFBP3-processing enzyme[40] might be related to growth characteristics. Consistent with retrogene analysis (Supplementary Note), there are multiple non-functional single-exon protease-like pseudogenes. Seven of these have complete ORFs without identified transcripts, indicating that they arose from recent retrotranscription events.

PRDM9, which binds DNA in recombination hot spots and affects recombination activity during meiosis[41] (Supplementary Note) is duplicated in catarrhine primates. Orthologs with all three functional PRDM9 domains have been computationally identified in placental mammals[42], however, these genes are often not in syntenic locations. In primates (including human and marmoset), panda, pig, and elephant, there is a PRDM9-like gene flanked by a conserved syntenic block including the genes URAH and GAS8. This gene, located near 16q telomere in human, is labeled PRDM7 in catarrhine primates but PRDM9 in marmoset and non-primates. Another gene (labeled PRDM9 in catarrhine primates) is located between the cadherin genes CDH12 and CDH10 at human 5p14[43]. This gene is present in chimpanzee, orangutan and rhesus macaque but is absent in marmoset and non-primates. The marmoset genome sequence provides two types of evidence that support the occurrence of a duplication in the catarrhine lineage after its divergence from platyrrhine primates: the phylogeny of PRDM9-like genes (Supplementary Figure 10b) and their genomic locations.

Population Genetics and Polymorphism

Genome sequence diversity was examined in nine marmosets (two animals from the New England Regional Primate Research Center (RPRC), two from Wisconsin National PRC (NPRC), and five from Southwest NPRC) (Supplementary Figure 11). This sample size is sufficient to identify common polymorphisms in this species, but will not detect a large proportion of low frequency or rare variants. Chimerism does not interfere with identifying SNPs that are polymorphic in the species as a whole, but does complicate the assignment of genotypes for specific SNPs to specific individuals. We investigated this effect by quantifying read balance (proportion of reads supporting each allele in apparent heterozygotes) and found different distributions in marmosets compared to a human control: more SNPs with read balance fractions between 5% and 25% are observed in marmosets. Simulations indicate that this flattened read balance distribution results from bases that are not polymorphic in the sampled animal but are either heterozygous or differently homozygous in the co-twin, with the low level of alternative reads representing the chimeric cells introduced during development (Supplementary Note, Supplementary Figure 2a).

We also explicitly modeled the expected numbers of sequencing reads covering a dimorphic SNP locus with one allele or the other, given a known fraction of chimerism and applied a maximum likelihood method to estimate the proportion of chimerism present in the marmoset samples from the sequencing data (see Supplementary Note). The chimerism fractions range from 12% to 37% (Supplementary Note, Table 6).

Using polymorphic autosomal bi-allelic SNPs (~7.7 million), we calculated pairwise allele-sharing genetic distances. To test if the genetic variation among individuals can be explained by their primate colony of origin, we performed principal component analysis (PCA) based on the pairwise distance. PCA separates the three colonies on the first 2 principal components (Supplementary Figure 11a) with individual M32784 from Southwest PRC more similar to individuals from other primate centers. Next, we used ADMIXTURE[44] to assess the ancestry of each individual. When K=3 (Supplementary Figure 11b), three groups corresponding to the colonies are identified. New England RPRC and Wisconsin NPRC individuals form distinct groups with little admixture. Consistent with the PCA result, two Southwest NPRC individuals (M32783 and M32784) show appreciable admixture from the other colonies (Supplementary Figure 11b). A neighbor-joining tree using the distance matrix (Supplementary Figure 11c) confirms the individuals from the same colony are grouped together, with the exception of M32784. The long terminal branch length suggests that most of the diversity exists among individuals.

We identified 107 polymorphic Alu insertions within common marmosets (Supplementary Figure 10a). Analysis of these insertions using Structure (version 3.3.2)[45,46] indicates population structure within the marmosets and detects two populations (Supplementary Table 41, Figure 12). The included marmosets show varying degrees of admixture with some individuals mostly assigned to one cluster and others assigned to both clusters (Supplementary Figure 12). The Structure analysis suggests the New England RPRC colony is assigned primarily to one cluster and the Wisconsin and Southwest NPRC colonies fall into the other cluster.

Discussion

Prior analyses of primate genomes have identified few specific changes that account for phenotypic differences among species, with the exception of genes that influence human brain size[47], language (reviewed[48]) or other unique human traits[49]. In contrast, our analysis presents a number of specific differences in gene content, miRNA number and sequence, and protein coding gene sequences in genes known to influence growth, reproduction and twinning propensity all potentially related to marmoset phenotypic adaptations (Supplementary Figure 13). Such divergence at multiple levels does indeed underscore the remarkable nature of this platyrrhine monkey species.

Online Methods

Genome sequencing and assembly

The 26.7M sequence reads, generated on ABI3730 instruments (Supplementary Table 1), average read length 700 bases (>=Phred[51] quality 20), were assembled using PCAP[52]. Assembly filtered for known non-marmoset sequence contaminants and singleton contigs and supercontigs < 2 kb included 99.98% of the input reads and is 59% A/T. WUGC Callithrix jacchus-3.2 was submitted to Genbank (UCSC version calJac3) and used by Ensembl to build gene models. Statistics (Supplementary Table 2) are for the initial assembly, before integrating finished BACs and adding interscaffold gaps and gaps representing centromeres and telomeres. The final assembly spans 2.91 Gb, with 2.77 Gb ordered and oriented along specific chromosomes (see below). The assembly represents an arbitrary consensus of the individual animal’s alleles. The non-repetitive assembly data were aligned against the repeat masked human genome at UCSC utilizing BLASTZ[39]. Orthologous and paralogous alignments[53] were differentiated and only “reciprocal best” alignments were retained and used to generate. the marmoset AGP files similar to methods described[12]. Documented inversions based on FISH data (see URLs and below) and inversions suggested by the assembly and supported by additional mapping data (e.g. fosmid and BAC-end pairs) were also introduced. Centromeres were placed based on their positions identified from cytogenetic data (see Supplementary Note). A total of 81 finished CHORI-259 marmoset BACs (totaling 15,576,643 bases) were merged into the final chromosomal files.

WUGSC sequenced Marmoset cDNAs (Supplementary Table 4) with Roche 454 Life Sciences instruments and methods[54] assembled using Newbler[55] and BLAT[56] aligned against the marmoset genome.

Using >700 human BAC clones we established synteny block organization of marmoset chromosomes and disambiguated inconsistencies and uncertainties of the genome assembly (see URLs to explore the data).

Gene Feature Annotation

The RefSeq[57] and Ensembl[18, 58] annotations used the general methods described (see URLs). Ensembl annotation raw compute stage (Supplementary Figures 6a) screened genomic sequence using RepeatMasker[59] (version 3.2.5 with parameters ‘-nolow -species homo –s’), Dust (Kuzio J, Tatusov R, and Lipman DJ: Unpublished but briefly described in[60])(together masking 47%) and TRF[61].

Predicted features include transcription start sites (Eponine–scan[62], FirstEF[63]), CpG islands (described[64]) and tRNAs[65]. Genscan results on RepeatMasked sequence were input for UniProt[66], UniGene[67] and Vertebrate RNA (see URLs) by WU-BLAST[68,69] alignments, resulting in 252,582 UniProt, 316,384 UniGene and 317,679 Vertebrate RNA sequences aligning.

Genewise[70] and Exonerate[71] produced coding sequence models using marmoset and human UniProt, SwissProt/TrEMBL (see URLs) and RefSeq[72] proteins mapped to the genome (Pmatch, Durbin, unpublished) (Supplementary Figures 6b, 6c). One model per locus was selected using BestTargeted module. Species-specific (here marmoset and human) data generated 1,908 (of 3,153) marmoset protein and 20,735 (of 22,320) human protein “Targeted stage” models with UTRs.

Raw Compute UniProt alignments were filtered, sequences with UniProt’s Protein Existence (PE) classification, level 1 and 2 were mapped with WU-BLAST and coding models built with Genewise in regions outside of Targeted Stage models, generating 57,019 mammalian and 42,323 non-mammalian additional “Similarity stage” models.

Marmoset cDNAs and ESTs and human cDNAs from ENA/Genbank/DDBJ with polyA tails removed were aligned to the genome using Exonerate[72] (Supplementary Figures 6d, 6e, 6f). With 90% coverage and 80% identity cut-offs, 139,713 (of 292,329) human cDNAs, 887 (of 986) marmoset cDNAs, and 2,562 (of 2,605) marmoset ESTs aligned. EST-based gene models (similar to those for human[73]) are displayed in a separate website track from the Ensembl gene set.

Similarity stage coding models were filtered to remove models with little cDNA or EST support, visualized using Apollo[74], extended using human cDNA and marmoset expressed sequences, resulting in 1,501 (of 2,119) marmoset, 13,150 (of 20,735) human, and 22,897 (of 31,863) UniProt coding models with UTR. Redundant transcript models were removed and remaining models clustered wherever any coding exons from two transcripts overlapped.

More information on the Ensembl automatic gene annotation process[19,20] is available in the references and the Supplementary Note.

Segmental duplications in Callithrix jacchus-3.2 were estimated using two computational methods; one compares assembly segments using BLAST (Whole Genome Assembly Comparison), WGAC,[8] the second, assessed excess depth-of-coverage of WGS data mapped to the assembly, WSSD[9]. All scaffolds were RepeatMasked (see URLs) and window masked[75] using the specific marmoset repeat library (Supplementary Note) composed of retrotransposons and other low complexity sequences. WGAC, identifies pairwise alignments >1 Kb and >90% identity. WSSD, identifies SDs > 10 Kb in length and >94% identity. For WSSD we mapped reads using Megablast with >94% sequence identity, >200bp non-repeatmasked bp and at least 200 bp of PhredQ>30 bp.

FISH Assessment of 97 clones on lymphoblast cell line nuclei and metaphases chromosomes from a marmoset unrelated to animal #186/17066. Duplicated probes had > 2 signals in 95-98% of > 60 observed nuclei (Supplementary Figure 3c). Sixteen clones showing strong hybridization background were tested three times without a clear pattern emerging and removed from further analysis. This unusual background may be due to incomplete masking by Repeatmasker and/or competitive hybridization conditions during FISH. Nine (of 16) belonged to the WGAC-WSSD+ category consistent with collapsed repeats.

As with ape assessment[76], We aligned 27,615,086 marmoset reads to human (build35, excluding random sequences) with repeat content masked (<20% divergent from consensus (RepeatMasker in either human or marmoset)). Aligned reads had >200bp of high quality (Phred score >27), >300bp aligned, >40% read length aligned, and <200bp repeat content. After evaluation, we applied an identity threshold of 85%, similar to the criteria applied in the macaque analysis. See Supplementary Note for details.

Sequence elements constrained in anthropoid primates were defined using the pipeline briefly outlined in Supplementary Note, described in detail[17]. To validate the functional role of the bioinformatically defined elements as transcriptional enhancers, we tested eight non-coding ASCs in embryonic stem cell (ESC) enhancer assays. Candidates were selected based on DNaseI hypersensitivity in human ESCs[77]. The eight human sequences and their mouse orthologs (identified using liftOver, Supplementary Table 11) were amplified from their respective genomic DNA, cloned into the SalI site, downstream of luciferase, of the pGL3-Pou5f1 vector using the Gateway Cloning System (Invitrogen), and transfected with the reporter constructs into hESCs (H1-WA-01, WiCell Research Institute, Inc.) and mESCs (E14TG2A, ATCC CRL-1821) using Fugene HD (Roche) or Lipofectamine 2000 (Invitrogen) respectively. Both cell lines are routinely tested for Mycoplasma contamination (Lonza Detection kit :LT07-318). A Renilla luciferase plasmid (pRL-SV40, Promega) was co-transfected as an internal control. Cells were harvested 48 hours after transfection and the luciferase activities of the cell lysate were measured using the Stop-Glow Dual Luciferase Reporter Assay System (Promega) (Supplementary Note).

MicroRNAs (877, Supplementary Table S2) were identified as being expressed or predicted based on cross-species conservation of mature microRNA or hairpin sequences.. Small RNAs were sequenced from total RNA from prefrontal cortex brain samples (A07-716monkB, 3.2 years, male; A09-122monkB, 12.8 years, female; A08-206monkB, 13.4 years, male; A08-337monkB, 13.0 years, female), and two placenta samples, using 36 bp reads on the Illumina 1G Genome Analyzer[78]. Usable reads were identified as described[78,79], omitting reads with <4 copies, <10 nt, >10 repetitive nucleotides, or that matched E. coli using WU BLAST[69] (Supplementary Table S2). Expressed microRNAs 100% conserved (Group A, 291 miRNAs) or with 1-3 mismatches (Group B, 240 miRNAs) in at least one other species in miRBase 17.0[80]. were identified. Known microRNAs in miRBase[80] 17.0 mapped to the marmoset genome identifed conserved (100% match, Group C, 119) or novel (with 1-3 mismatches, Group D, 120) miRNAs.. Group A to D sequences (~22nt) aligned with BLAT (−stepSize=5 repMatch=100000 - minScore=0 −minIdentity=0 −fine), and their flanking sequences (+/− 200 bp) extracted from UCSC were folded using Vienna RNAfold[78] twice to confirm hairpin structures with the mapped sequenced in the mature miRNA location. Group E contains the 91 novel miRNAs (20 passed high stringency filters) identified, trimmed to include only the hairpin bases (60–150 nt) (Supplementary Table S2).

WU BLAST comparison identified marmoset micoRNAs conserved in four anthropoid primates (−nogaps −N −1000 −mformat=2 −warning −kap −hspmax=10) (marmoset: calJac3, human: hg18, rhesus: rhemac2, orangutan: ponAbe2, chimpanzee: panTro2 from UCSC). BLAT mapping (−stepSize=5 −repMatch=100000 −minScore=0 −minIdentity=0 –fine) of precursor miRNA hairpins on marmoset Chr. 22 to rhesus, orangutan and chimpanzee identified the best matches which were realigned to marmoset miRNAs hairpins using Smith-Waterman to identify nucleotide changes in the mature microRNA sequence. Human Chr. 19 hairpins were mapped to calJac 3 using Galaxy lift over and BLAT alignment and realigned as above (see conservation in Supplementary Tables S3, S4, S5, S6, S7, and S8).

MicroRNAs predicted using SVM (Group F)

Human precursor miRNAs (miRBase 14.0[81]) WU-BLASTN[68,69] (see URLs) matches > 20 bp to calJac 3.2. (−M 1 −N −1 −Q 3 −R 2 −W 9 −filter dust −mformat 2 −hspsepSmax 40 −e 1e-3) were extended to match their entire length and realigned using MAFFT[82] (maxiterate 1000 --localpair –quiet). Matches were identified with (i) length > 40 bp, (ii) 100% conserved seed region (mature miRNA nucleotides 2-8), (iii) > 90% mature miRNA sequence identity, (iv) total precursor conservation > 50% length, (v) at most 2 gaps in mature miRNA, (vi) minimum free folding energy (MFE) < −15 kcal/mol, (vii) > 40% of bases paired, (viii) mature regions not overlapping a multiple loop region, and (ix) probability < 5% for a randomly shuffled hit sequence to have a lower MFE than the native sequences for <95% conserved matches. The lowest e-value hit for overlapping loci was subjected to a Support Vector Machine (SVM) model trained to distinguish miRNAs from unspecific genomic stem-loop sequences or other non-coding RNAs. Developed for miROrtho annotation database[83] (see URLs), the model incorporates themodynamical, structural and sequence features found in known miRNA genes. Using an initial BLAST e-value cutoff of 1e-6, SVM score over 0.5 and 100% mature miRNA sequence conservation to any known miRBase miRNA identified 589 genes (Group F)..

small RNA sequences (above) mapping within 4 bp on the same chromosome as the microRNA, normalized by total usable reads. Euclidean hierarchical clustering of genes and arrays with Cluster 3.0 and Tree View[84] (see URLs) used the log2 transformation of microRNAs per 10 million usable reads with the median expression value across the six samples set to zero.

MiRmap[85] identified mRNAs with 3′ UTR matches to miRNA bases 2 - 8 and predicted the repression strength with a model encompassing thermodynamic, conservation, probabilistic, and sequence-based approaches. We computed the total energy of the miRNA-mRNA duplex (similar to[86]), branch length score[87] and SPH test implemented in PhyloP[88] and computed statistical significance of the seed-match based on 3′-UTR sequence composition. The three features of TargetScan context score[89] were included in miRmap for a total of 11 features, of which 3 are novel (see URLs). These data were generated by mapping all human RefSeq genes to marmoset based on UCSC “Other RefSeq” track, multiple mapping locations in marmoset were retained, represented {refseqAccession}.1, {refseqAccession}.2, etc. Where the 3′ UTR differs between different mapped locations, this could represent true paralogs or assembly errors. The extracted marmoset 3′-UTRs were aligned using MAFFT[82] to the TargetScan 5.1 23 way UTR alignments and marmoset target genes identified with 3′-UTR binding sites to the mature marmoset Chr. 22 family miRNAs.

Conservative 1:1 orthologs between marmoset and human, chimpanzee, rhesus macaque, orangutan, mouse, rat, and dog produced using UCSC[90] whole-genome alignments and genes (July 2010) including partial transcripts missing 10% on both ends. Transcripts on chromosomes > 100 nucleotides from RefSeq (58,126) , knownGene (118,345), ENSEMBL (128,193) , and VEGA (73,873) clustered into 21,694 genes based on location.

Each transcript was transferred to other species and subjected to tests designed to exclude genes that have undergone large-scale changes other than point mutations (as in[19]), testing for breaks in synteny, significant assembly gaps overlapping the transcript, frameshift and non-sense mutations, conservation of gene structure elements (splice sites, start codons, stop codons), and recent duplications causing mis-assignment of 1:1 orthology. Clean transcripts pass all tests. We choose a representative clean transcript for each locus, preferring longer transcripts clean in more species (summarized in Supplementary Table 12). This conservative set (13,717 1:1 orthologs of human and marmoset), has 41% covering all eight species, 27% missing one species, 15% missing two species, 10% missing three species, and less than 7% missing more than three species.

Gene family evolution was investigated in four other primates, two rodents, and three Laurasiatheria with fully sequenced genomes (human, chimpanzee, orangutan, rhesus macaque, marmoset, mouse, rat, dog, horse and cow). Gene families, including gene and protein names and genome coordinates were retrieved from Ensembl gene trees, version 58 (see URLs). Genes with multiple short introns (< 50 bp), short coding regions (< 100 bp), or present in < 3 species were removed, and we analyzed separately families with genes in only one lineage (Euarchonta, Glires, Laurasiatheria). The final set included most genes and families from the original Ensembl annotations (Supplementary Table 13), and was used to infer the ancestral family size with maximum-likelihood CAFE[91] analysis with the following ultrametric tree built according to[92]: ((((((chimp:6,human:6):7, orang:13):11, macaca:24):16, marmoset:40):47, (mouse:17,rat:17):70):6, ((dog:74,horse:74):9,cow:83):10). where numbers correspond to million years (Supplementary Note).

Positively selected genes in the 1:1 orthologs were identified using Markov models of codon evolution and maximum likelihood methods similar to PAML[93]. Further downstream analysis such as enrichment analysis for GO categories was performed as described[19]. Supplementary Note details genes identified using FDR<0.01.

Genes involved in growth pathways and twinning

Candidate genes identified using 33-way EPO alignments[18] containing marmoset NS (compared to human) conserved in haplorhine primates (human, chimpanzee, gorilla, orangutan, rhesus macaque and tarsier) were sequenced. The NS effect was defined using SIFT[94] and some candidates were omitted due to conflicting evidence. Genes and coordinates are listed in Supplementary Table 39. Species include: Saguinus bicolor martinsi*, Saguinus imperator imperator, *Saguinus midas niger, Saguinus fuscicollis weddelli, Callithrix cebuella pygmaea*, Leontopithecus rosalia*, Cebus apella, Callimico goeldii, Ateles belzebuth, and Saimiri sciureus (species with an asterisk were also selected for miRNA sequencing). Sanger sequences were assembled (Velvet[95]), mapped to the genome (BLAT[51]), and aligned (MAFFT[82]). In 49 of 82 exons sequenced, data was insufficient to determine if the marmoset NS were callitrichine or NWM specific (Supplementary Note).

Protease Genes

We mined the marmoset genome for protease genes (see URLs) using BATI (Blast, Annotate, Tune, Iterate). Curated human proteases were compared to the marmoset genome with TBLASTN algorithm using the tbex script and the locations of marmoset protease genes predicted with bsniffer. Putative novel proteases were predicted with bgmix (Supplementary Note) and visually inspected.

Variation Analysis

SNPs (7,697,538) in reads aligned to the genome using BWA (version 0.5.9-r16, default parameters), were called using SAMTOOLS[96] (version 0.1.14 (r933:176), command “$ samtools pileup −Bvcf $ref_genome $bam”), filtered q>20, D<100), with monomorphic, multi-allelic, and singletons removed. Pairwise allelesharing genetic distance was calculated[97] and the resulting matrix used for PCA and neighbor-joining tree construction (MATLAB ver. r2010b). Genetic ancestries of each individual was determined with ADMIXTURE[44] in a given number of populations without using population designation. We filtered out SNPs with r2 > 0.2 within 100 SNP window using PLINK[98], leaving 411,924 autosomal SNPs.

Alu genetic analysis

Best matching loci from CalJac3.2 for each Alu subfamily were identified using BLAT[51] or retrieved from a local RepeatMasker (see URLs) analysis using a custom library and subfamilies with evidence of recent mobilization (divergence up to 1%) from the consensus sequence used for population genetic analyses. For phylogenetic analyses, Alu insertions of subfamilies were selected with varying divergence from the consensus sequence.

We retrieved marmoset Alu elements with ~500 bp flanking sequence, identified orthologous loci using BLAT[51] and retrieved the sequences if the flanking sequence matched unambiguously in the other genome and the Alu insertion was absent. We did this for human, chimpanzee, orangutan, and rhesus macaque. We aligned the flanking sequence (BioLign/BioEdit) and selected primers (Primer3[99] or selected manually) to minimize the nucleotide substitutions and other Alu insertions. Primers were tested using UCSC In-Silico PCR[51] and synthesized by Sigma Aldrich.

PCR amplifications (96-well format) in 25 μl using Perkin Elmer GeneAmp 9700 or a BioRad i-cycler thermo-cycler contained 15-25 ng of template DNA; 200 nM of each primer; 1.5-2 mM MgCl2, 1X PCR buffer (50 mM KCl; 10 mM TrisHCl, pH 8.3); 0.2 mM dNTPs; and 1-2 U Taq DNA polymerase. Conditions included: initial denaturation at 94 °C for 90 sec, followed by 32 cycles of denaturation at 94 °C for 20 sec, 57 °C annealing for 20 sec (see URLs for exceptions), and 72 °C extension for 30 to 70 sec depending on amplicon size, with a final 72 °C extension for 2 min. If necessary, temperature gradient using HeLa DNA identified the optimal annealing temperature. 20 μl of each reaction was fractionated in 2% agarose gel containing 0.1 μg/ml ethidium bromide at 175V for 50-60 minutes and visualized with UV-fluorescence.

Using genotype data from unlinked markers we inferred the population structure omitting information of the origin of the samples with a model-based clustering analysis[45,46] under the admixture model that assumes that individuals may have mixed ancestry.

The number of identifiable population clusters (K) with the highest likelihood was determined using initial values of K (1 to 5), burn-in period (1,000,000 iterations) and run-length of 1,000,000 steps repeated at least five times. After determination of K (2) 25 replications were run under identical burn-in and run-length settings. The Structure analyses were run on a desktop machine with 4 CPUs.

Animal samples used in these study were obtained under protocols approved by the relevant institutional animal care and use committees from animals maintained in AAALAC-accredited animal care programs.