The pioneering role of PRDM9 indel mutations in tarsier evolution

PRDM9 is currently the sole speciation gene found in vertebrates causing hybrid sterility probably due to incompatible alleles. Its role in defining the double strand break loci during the meiotic prophase I is crucial for proper chromosome segregation. Therefore, the rapid turnover of the loci determining zinc finger array seems to be causative for incompatibilities. We here investigated the zinc finger domain-containing exon of PRDM9 in 23 tarsiers. Tarsiers, the most basal extant haplorhine primates, exhibit two frameshifting indels at the 5′-end of the array. The first mutation event interrupts the reading frame and function while the second compensates both. The fixation of this allele variant in tarsiers led to hypothesize that de- and reactivation of the zinc finger domain drove the speciation in early haplorhine or tarsiiform primates. Moreover, the high allelic diversity within Tarsius points to multiple effects of genetic drift reflecting their phylogeographic history since the Miocene.

zinc fingers of anthropoid primates 7 and tarsiers, and discussed the possible function of PRDM9 as driving force behind the anthropoid-tarsiiform split and the divergence of the tarsiiform lineage.

Results
Tarsier PRDM9 sequence characterization. The examined exon of Western, Philippine and Sulawesi tarsiers was homologous to exon 11 of the human PRDM9 gene and was about 1037-1793 bp long depending on the number of zinc finger repeats. It could be divided into two parts, the 5′ sequence and the 3′ C2H2 zinc finger array. Structural and functional properties of the zinc finger protein encoding exon are shown in Fig. 1.
In contrast to the variable C2H2 zinc finger domain which was used for the analyses, the 5′ sequence of the exon was more conserved. This part of the sequence contained a C2H2 zinc finger (ZnF1) at the beginning of the exon and four presumably former functional zinc fingers. The second and third zinc fingers were both degenerated by the loss of a zinc ion binding cysteine ligand, see Fig. 2, which is in principal not uncommon in zinc finger  Shown are aligned nucleic acid and peptide sequences comprising zinc fingers 3-6 (ZnF3-ZnF6) of Homo sapiens (HSA) and Tarsius (represented by a Sulawesian individual, CD26) and illustrating the frameshifting nature of insertion (blue) and deletion (green) mutations. Note the nonsynonymous substitution in the very first triplet which turns cysteine to serine (orange) and degrades ZnF3 of the tarsier. Amino acid sequences of human ZnF3-ZnF6 and ZnF6 of the tarsier show the classical C2H2 zinc finger motif (indicated by bold and underlined amino acids).
Scientific RepoRts | 6:34618 | DOI: 10.1038/srep34618 arrays 22 . The sequence of the third zinc finger was also interrupted by a 2 bp-insertion that caused a framing error affecting the next two zinc fingers which we refer to as former zinc fingers. A deletion in the second altered zinc finger, which would be the fifth overall, restored the reading frame. The array of classical 84 bp C2H2 zinc fingers started with the end of this fifth zinc finger (Figs 1 and 2).
Indel events are quite rare in coding regions due to purifying selection 23,24 , especially if their length is not a multiple of three 25 . While the insertion results in a different amino acid sequence and may alter protein function, the deletion alone generates premature stop codons which lead to a faster pseudogenization respectively gene loss 26 . We therefore consider it likely that the insertion predates the deletion mutation which then in turn acted as a compensatory mutation. Compensatory mutations are twice as common as reversions to the original state and are often found in close proximity of the initial mutation 27,28 being consistent with our findings. As only both indels jointly restore protein integrity we further assume that the temporal offset between indel mutations resulted in a temporary loss or limited protein function.
Allelic variation of the zinc finger array in tarsiers. The alleles are defined on amino acid level, i.e. synonymous nucleotide changes are excluded, and are restricted to the C2H2 array. We found 28 alleles in 23 individuals with 15 individuals being heterozygous, (Table 1, Fig. 3).
The 25 alleles exclusively found in Sulawesi tarsiers show similarities, most often within populations (see Materials and Methods for sample sites). They mainly differ in the number of zinc fingers but hardly in their sequence like the alleles 9 and 10. Alleles 4, 5 and 6 have each ten zinc fingers with minor changes in the amino acid sequence (Fig. 3). Only allele 14 is shared between two Sulawesian populations (Luwuk and Laone), both belonging to Tarsius dentatus.
There are two types of 5′ -most C2H2 zinc fingers differing in one binding amino acid (HHR, HRR). One variant can be found on the northern peninsula (Duasaudara, Labanu, Ogatemuku) and in Korosule. The other is specific to populations inhabiting the central, eastern and southern parts of Sulawesi, with the exception of one copy also occurring in north-eastern Sulawesi. The 5′ -most C2H2 zinc fingers accord to some extend with the two lineages found in Driller et al. 21 . Lineage 1 comprises the populations on the northern peninsula and Kendari, lineage 2 all other populations.
The 3′ -most C2H2 zinc fingers are, with one exception, identical regarding the three binding amino acids RHT (Fig. 3). Five of them show a nonsynonymous substitution which, however, does not affect the DNA-binding sites. They were observed in the populations of Labanu, Korosule and Koja. Three out of 19 zinc fingers with key codon positions QSR, LSR, and RNT are unique. Some are only found in few populations like the motifs specifying the binding triplet QNI or RNR and being restricted to Laone and Kendari, or Labanu and Ogatemuku, respectively. Especially the hints on the relationship of the populations on the northern peninsula among each other and between Laone and Luwuk are comparable with results of Driller et al. 21 . The relatively conserved 5′ -most and 3′ -most C2H2 zinc finger, the more variable ZnFs in between and the appearance of similar ZnFs in different species are also seen in mice 29 .
Alleles of T. bancanus and T. syrichta are species-specific encoding zinc finger motifs not present in Sulawesi tarsiers. Despite their long geographic isolation they still share two zinc fingers (Fig. 3). Both alleles of the Western tarsier end with a zinc finger degenerated by the exchange of the second zinc binding histidine ligand with asparagine. With only three repeats the zinc finger domain of the homozygous Philippine specimen is unusually short, but according to Stubbs et al. 30 capable to bind DNA. However, more individuals or rather more alleles of both species could reveal further and better insights regarding the relationship of PRDM9 zinc fingers between the three main tarsier groups.

Phylogenic analysis.
To estimate the phylogenetic affiliations of the sequenced tarsier zinc fingers within the haplorhine divergence we reconstructed a phylogeny with Microcebus murinus as strepsirhine outgroup. Including all first degenerated and C2H2 zinc fingers (Fig. 4a) it completes the results of Schwartz et al. 7 from a haplorhine perspective. The basal trifurcation separates degenerated zinc fingers of tarsiers, degenerated zinc fingers of anthropoid primates and all C2H2 zinc fingers with node support values above 0.9 (Fig. 4a). Within the C2H2 cluster tarsier zinc fingers form a well-supported monophyletic group. Within this group 5′ most and 3′ most C2H2 zinc fingers are clustered reflecting the biogeography of the individuals tested. Apart from these findings our tree topologies (Fig. 4) are congruent with those obtained for anthropoids as presented by Schwartz et al. 7 .
The second phylogeny (Fig. 4b) is based exclusively on first degenerated zinc fingers. Here again tarsiers constitute a monophyletic group basal to the anthropoid clade. Within Tarsius Sulawesi tarsiers are monophyletic with respect to Western and Philippine tarsiers.
Positive and negative selected sites. Zinc finger genes are known to alter quickly in their zinc finger array and in particular at the DNA-binding codons (positions − 1, 3 and 6 relative to the α -helix) which show signatures of positive selection [31][32][33] . The PRDM9 zinc finger domain is no exception as shown by positively selected DNA-binding amino acids in rodents and primates 3,34,35 . The high mutation rates at codon positions directing the interaction of PRDM9 with DNA is deemed responsible for the formation of new hotspots and thus antagonizing hotspot erosion 5,36,37 .  Table 1. Each zinc finger is depicted with a triplet composed of its binding amino acids (− 1, 3, 6). Each amino acid is marked with a distinct colour. Dots hint towards one or more nonsynonymous substitutions apart from the shown amino acids. Abbreviations: T.s., Tarsius syrichta; T.b., Tarsius  We found positive selection at the three DNA-binding codons, i.e. positions − 1, 3, 6 relative to the α -helix (see Fig. 5). Structural important sites like the zinc ion ligands are conserved by negative selection. We consider the ongoing selective pressure -positive as well as negative -acting on tarsier PRDM9 as relevant proof for gene integrity being important in the context of a probable temporary loss of function which was subsequently abolished by a compensatory mutation.

Discussion
With four non-functional zinc fingers at the 5′ end of the tarsier zinc finger array (Fig. 1), degenerated by missing cysteine ligands and frameshifting indels, we detected a completely new PRDM9 allelic variant in primates and mammals in general. This sequence singularity is also reflected in our inferred phylogenies including first degenerated and all functional C2H2 zinc fingers of the array and clearly designating tarsiers as monophyletic group within haplorhine primates (Fig. 4a,b). All these sequence autapomorphies, but especially those associated with a temporary loss of allelic integrity (see results) might possibly indicate an active role of PRDM9 in the divergence process between anthropoid and tarsiiform primates or alternatively, along the tarsiid lineage. A possible speciation scenario could be as follows: In a population of haplorhine or tarsiid ancestors a damaging, frameshifting mutation occurs in one allele (Figs 2 and 6). Heterozygous individuals may be subfertile because gene function is reduced due to one non-functional allele. However, Brick et al. 6 showed that one allele can determine 75% of the hotspots while the other is responsible for 22% [1% of all hotspots are shared and 2% are new in comparison to parental homozygote constellations]. Furthermore, male mice which are heterozygote null for Prdm9 are not sterile but reproduce at a later stage for the first time and father fewer offspring, both in comparison to the wild-type 38 . Individuals homozygous for the damaged allele are probably infertile due to a considerably reduced or loss of gene function and comparable to Prdm9 knockout mice where a meiotic arrest at pachytene stage leads to sterility 4 . The second frameshifting mutation restores the functionality of the former erroneous allele (Figs 2 and 6) and is therefore positively selected and fixated 39,40 .
Two of the now three possible allelic variants in homozygotes are functional (Fig. 6), the one with two original or wild-type alleles and the other with two alleles carrying both indel mutations. By contrast, all heterozygotes shown in Fig. 6 are possibly less fit or rather subfertile resulting from a decreased recombination activity affected by either one damaged allele or different hotspot usage among the two types of functional PRDM9 alleles 41 . The novel allele finally get fixed in a precursor of tarsiiform primates but in any case no later than 20 MYA, the time to the most recent common ancestor (MRCA) of the three extant tarsier clades 17,19,21 , all exclusively carrying the new allelic variant. Explanations for the fixation of the mutant allele in an ancestral population of modern tarsiers may include directed and/or random evolutionary processes like biased gene conversion or genetic drift. With regard to the former it has been shown that a newly introduced PRDM9 allele increased hotspot activity in mice through altered and haplotype-biased DNA-binding, thus counteracting fitness loss during the process of hotspot erosion by allelic gene conversion 36 . The mutant allele could have had a similar hotspot enhancing and therewith beneficial fitness effect in a population of ancient tarsiids or their progenitors. Given that its selective advantage was sufficiently large the novel allele could also have overcome the effects of genetic drift where typically rare beneficial mutations are more likely to go extinct 42,43 . On the other hand, we could suppose that it was this process of random change in genetic composition that maintained the newly arisen allele resulting in the structural uniformity of the PRDM9 zinc finger domain and facilitating divergence of an ancestral lineage finally leading to modern tarsiers.
Extinct relatives of living tarsiids branched away from anthropoid primates very shortly after the haplorhinestrepsirhine split 9,19,44 . Depending on the data, either paleontological or molecular, first primates evolved sometime between 55-87 MYA 9,15,19,45 . Within this period two drastic climate changes occurred and both induced mass extinctions 44,46,47 creating conditions promoting adaptive radiation 44 . Proto-tarsiiforms could have survived or emerged from these evolutionary bottlenecks possibly also because PRDM9 produced new recombination landscapes allowing allele combinations more favourable or adaptive to the changing environmental conditions and/ or newly vacant ecological niches. Another extinction event associated with Eocene-Oligocene cooling 48 largely decreased primate diversity with the disappearance of omomyids, a controversial fossil tarsiiform primate [49][50][51] . This event presumably caused a geographic shift in tarsier distribution from mainland to insular Southeast Asia 52,53 exposing tarsiids of modern aspects to another population bottleneck and thus providing a further option for the fixation of the novel tarsier-specific PRDM9 allele. A heightened relevance of genetic drift in tarsier molecular evolution was already discussed 54 , as extant tarsiers, and especially those endemic to the Indonesian island of Sulawesi, have been subject to significant climatic and tectonic change 55,56 that triggered allopatric speciation in the Malay Archipelago since the Miocene 17,18,20,21,[56][57][58] . Hence, the much longer history of independent evolution in environmental instability may gives more credence to the theory that the salient molecular changes of the PRDM9 zinc finger array described here have developed along the tarsiid lineage rather than initialized the anthropoid-tarsiid split.
PRDM9 variation in alleles, zinc finger motifs and especially at the key codon positions reflect the different levels of evolutionary independence of the three major tarsier clades (Western, Philippine, and Sulawesi tarsiers), provided the single non-Sulawesi specimens are representative for the respective biogeographic region. Like it has been observed in mice, chimpanzee and humans, 5′ -most and 3′ -most C2H2 zinc fingers are more conserved or even identical within or at least between closely related species 29,41,59,60 even if positive selected sites are factored in. In tarsiers this appears to apply mainly to species endemic to the same biogeographic region, but is also fairly reflected in lineage-specific constraints on 5′ -most C2H2 zinc fingers of Sulawesi tarsiers 21 . Considering the three DNA binding amino acids at internal zinc fingers of Sulawesi tarsiers we found several duplicated zinc finger motifs in each species/subspecies indicating their recent origin 61 for most populations being estimated at less than 500,000 years ago 21 . Particularly striking are the blocks of zinc finger motifs in northern Sulawesi tarsier populations (Fig. 3). The order of their occurrence very plastically depicts how these populations very recently evolved from a single common ancestor 21 .
In conclusion, the high mutation rate of the PRDM9 zinc finger domain 3,7,29 and multiple events of genetic drift produced an enormous diversity of PRDM9 alleles in tarsiers. The geographic distribution of zinc finger alleles and especially the occurrence and enrichment of specific zinc finger motifs reflect phylogeographic patterns of extant tarsiers further strengthening an involvement of PRDM9 in population differentiation. More intriguing, however, is our discovery of a hitherto unknown and, in addition, tarsier-specific PRDM9 allele variant that arose from two indel mutations and played a perhaps decisive role in tarsiid evolution.

Material and Methods
Sample set. The sampling of 23 individuals comprises 21 specimens of nine Sulawesi tarsier populations (1-4 individuals/population, see Fig. 7) and two non-Sulawesi tarsiers (Tarsius bancanus and T. syrichta). Samples from Sulawesi were obtained from previous studies 18,21 , while the Western and the Philippine specimen were provided by Y. Rumpler (Les Hôpitaux Universitaires de Strasbourg, France) and J. Brosius (University of Muenster, Germany), respectively. Whole genome amplifications (WGA) of each sample (40-80 ng/μ l) were used as a template in PCR to amplify the exonic region of the PRDM9 gene containing the zinc finger domain. PCR and Sequencing. We PCR-amplified and sequenced the exon which encodes the zinc finger domain of the PRDM9 gene in three parts: the 5′ and 3′ flanking regions and the repetitive elements in between (see Table 2).
To reduce amplification of unspecific products we conducted wax-mediated hotstart PCR. Each PCR-reaction contained 30 μ l with final concentrations of 200 μ M dNTPs, 2.5 units Taq DNA polymerase, 1x PCR buffer including 1.5 mM MgCl 2 (QIAGEN Taq PCR Core Kit), and 0.33 pM per primer. PCRs were run under the following conditions: 3 minutes of initial template denaturation at 94 °C was followed by 35 cycles of denaturation (40 sec at 94 °C), primer annealing (1 min at primer-specific temperatures) and DNA elongation (1-1.5 min at 72 °C). A final elongation step of 5 minutes at 72 °C finished the PCR. Product sizes were estimated on ethidium bromide stained 1.5% agarose gels together with O'RangeRuler 100 bp DNA Ladder and GeneRuler 100 bp Plus DNA Ladder (both Fermentas/Thermo Scientific).
In general PCR products were enzyme purified before sequencing. In some individuals PCR reactions yielded multiple different sized products, either because of size variation between the two alleles or unspecific amplicons. The separation of the two alleles or the isolation of the product from non-specific sequences was done by gel extraction with the QIAquick Gel Extraction Kit (QIAGEN). Purified PCR products were sequenced on both strands using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) and the corresponding PCR primer pair. After an SDS/heat treatment for elimination of unincorporated dye terminators sequences were processed on an ABI 3130xl genetic analyzer (Applied Biosystems).
For individuals with ambiguous sites present in their DNA sequences PCR products were purified by ethanol precipitation, ligated into a plasmid vector (pGem-T Vector System I, Promega) and transformed into One Shot TOP10 Chemically Competent E. coli cells (Invitrogen). Where gel extraction failed to separate alleles of different lengths sequences were also isolated by cloning.  PRDM9-R2 CCTCATCTCTAGTCATGAAAGTGG [24] PRDM9-R3 GAAGCACCCTCCAAGCTG [18] PRDM9-R4 GGTCTCTTTACACTCTTGSAGTC [23]  At least six positive clones, three of each allele, were PCR amplified and sequenced. Colony PCR was carried out in 20 μ l reaction volumes according to standard protocols using the Taq PCR Core Kit (QIAGEN) and vector-specific primers. For long amplicons where sequencing read length was not sufficient we used the internal primer PRDM9-F2 (see Table 2) to get the middle portion of the target sequence. Sequence Analyses. Raw sequences were edited and consensus sequences (KU948327 -KU948367) were generated in BioEdit 7.1.3.0 62 . To reveal the structure and integrity of the tarsier PRDM9 gene we created alignments in BioEdit with the following approach. We compared 5′ -terminal sequences including the degenerated zinc fingers in front of the C2H2 array 1) within Tarsius, and 2) between tarsiers, anthropoids (Homo sapiens, ENST00000296682; Pan troglodytes, GU166820.1) and rodents (Rattus norvegicus, ENSRNOT00000066370; Mus musculus, ENSMUST00000167994). We further examined the zinc finger domain by comparison of isolated zinc fingers within and between taxa of Tarsius, anthropoids and rodents. Structural and functional properties of tarsier PRDM9, in particular with regard to the type of zinc finger motifs, was also validated by peptide sequences obtained with EMBOSS Transeq v. 6.3.1 63,64 . Zinc finger Phylogeny. We generated two tree topologies, one including all 5′ -most degenerated zinc fingers and one including both, functional C2H2-zinc fingers and degenerated first zinc fingers of the array. Haplorhine taxa represented in each data set comprised Tarsius (this study), Homo sapiens (ENST00000296682) and several other anthropoid primates, whose sequence data were adopted from Schwartz et al. 7 . Microcebus murinus served as outgroup. Corresponding sequence information was extracted from the Ensembl data base (assembly micMur1, database version 77) and by using BLAST+ 2.2.28 65 . The zinc fingers were isolated with the first codon for cysteine as starting point. All except 5′ -most degenerated and several 3′ -most zinc fingers had a length of 84 bp. Lacking the first cysteine-triplet 5′ -most degenerated zinc fingers start with the second codon for cysteine, and therefore being nine bp shorter. 3′ -most zinc fingers can lack base pairs when interrupted by a stop codon, but they contain at least both cysteine and histidine ligands. Two 3′ degenerated zinc fingers of Tarsius bancanus were included. Like Schwartz et al. 7 we excluded the four binding triplets. To reduce the calculation time we removed identical zinc fingers within genera with the Fabox DNAcollapser 66 . Models of nucleotide substitution were determined with jModelTest 2.1.6 67,68 . Phylogenetic trees were generated with the Maximum Likelihood method (PhyML 3.0 67 ) including an approximate likelihood ratio test (aLRT) and SH-like supports. The degenerated zinc finger of Microcebus murinus served as outgroup in both tree topologies. Nodes with support values below 0.5 were collapsed in TreeGraph 2.2.0-407 beta 69  Detecting selective pressure. All functional C2H2 zinc fingers (84 bp) were tested for selective pressure on codon sites. We inferred a phylogenetic tree based on all non-identical sequences determined by Fabox DNAcollapser 66 using the Maximum Likelihood method [PhyML 3.0, settings other than default: K80, bootstrap: 100] 67 . The best fitting substitution model was selected using jModelTest 2.1.6 67,68 . Tests of selection were performed with the sitewise-likelihood ratio method (SLR 71 ) using default settings.