Exploring the genetic mechanisms driving KIR diversification

Abstract

Killer cell immunoglobulin-like receptors (KIRs) are key modulators of natural killer cell activity, displaying either activating or inhibitory potential upon recognition of major histocompatibility complex (MHC) class I molecules. The genomic organization of KIR genes is complex, involving copy number variation and allelic polymorphism, which is probably due to their coevolution with highly polymorphic MHC ligands. The KIR diversity is reflected by more than 70 similar region configurations encountered in humans, generated through meiotic recombination events. Rhesus macaques happen to display even more diversity, and over 100 distinct configurations were identified in a relatively small cohort of animals. More than half of these region configurations feature hybrid KIR genes, suggesting a more pronounced mode of diversification in macaques. The molecular mechanism facilitating meiotic rearrangements in the KIR region is poorly understood. Examination of 21 rhesus macaque and 14 human KIR region configurations revealed the presence of long terminal repeats and PRDM9 binding motifs associated with recombination hotspots. The variable DNA recognition patterns of PRDM9 could potentially contribute to the differing recombination activities documented for the KIR region in humans and macaques. The diversification process of the KIR repertoire in natural killer cells is fundamentally distinct from the mechanisms generating T and B cell receptor diversity or MHC polymorphisms. This sophisticated recombination machinery preserves the functional integrity by the frequent generation of in-frame KIR genes. A diverse KIR repertoire contributes to the protection of individuals and populations against pathogen evasion and subversion.

Introduction

Natural killer (NK) cells are lymphocytes that contribute to the early response during infection or tumor formation. In contrast to T and B lymphocytes involved in the adaptive immune response, the innate counterparts do not feature rearranged antigen-presenting receptors but monitor their environment via an array of germline encoded molecules.^1–3 Important modulators of NK cell activity are the killer cell immunoglobulin-like receptors (KIRs), which exhibit activating or inhibitory potential upon recognition of polymorphic epitopes on major histocompatibility complex (MHC) class I molecules.⁴ Coevolution of these 2 immune families, driven by a constantly changing landscape of pathogens, propels their genetic conservation or diversification.^5,6 In total, 17 KIRs are documented in humans, which display diversity in receptor structure, signaling potential, ligand binding specificity and affinity, expression levels, and cellular localization.^7,8 This variation is extended by polymorphisms, as evidenced by a repertoire of 1,600 KIR allotypes.⁹ Only a limited set of KIRs are encoded within an individual’s genome, whereas the whole scale of allelic variation becomes visible at the population level. This strategy of the immune system likely minimizes the risk of an entire population being highly susceptible to a specific pathogen and reduces the potential for pathogen evasion. Even more, NK cells express 1 to 3 of these available KIRs in a stochastic manner, a process maintained by the epigenetic status of promoters, generating NK cell heterogeneity at an individual level.^10–13

At a genomic level, the plasticity of the KIR system is further augmented by extensive gene copy number variation. Over 70 different gene configurations (unique sets of genes present on 1 chromosome) are defined for the human KIR gene cluster on chromosome 19q13, the majority of which consist of a centromeric and a telomeric region.^14–16 The 7 most frequent KIR gene configurations are different combinations of 4 centromeric and 2 telomeric regions. The less common gene cluster organizations originate from chromosomal rearrangements that propel structural diversity by the introduction and deletion of genes, by recombining complete centromeric and telomeric regions, and by the reshuffling of gene segments. This latter molecular process may generate novel hybrid entities that might encode KIRs with distinct functional properties.^15,17,18 These recombination events make the KIR cluster one of the most variable regions in the human genome, comparable to the extensive genetic diversity recorded for their MHC ligands.

The rapid evolution of this immunogenetic system is also substantiated by the existence of diverse KIR region configurations and gene repertoires in other primate species.^19–22 For example, unparalleled KIR gene diversity is documented for rhesus macaques, a species that is relevant as model for infectious diseases in biomedical research.^18,19,23,24 Over 118 KIR haplotypes (gene configurations with allelic polymorphism taken into account as well), representing 102 region configurations, were recorded in a relatively small outbred colony of 300 rhesus macaques, indicating a highly unique KIR repertoire per individual.²⁴ The large number of distinct gene configurations most likely originates from frequent recombination events, a phenomenon that appears to be more pronounced in the KIR region of rhesus macaques as compared with its human counterpart. This difference might be explained by a stronger diversifying selection acting on the macaque KIR region, potentially driven by several factors, such as the imperative to maintain immune recognition in the face of evolving pathogens, reproductive success, and coevolution with the dynamic nature of MHC class I genes, which, in contrast to humans, display substantial levels of copy number variation in macaques as well.²⁵

The differential diversifying activity in the human and macaque KIR regions likely involve molecular aspects that promote recombination, but these putative features are poorly understood. Across the genome most recombination events are confined to narrow hotspots that are unevenly distributed on chromosomes.²⁶ A higher hotspot density is in general observed at the terminal subtelomeric regions,²⁷ where indeed on chromosome 19 the KIR gene cluster is situated in primate species. In addition to chromosomal location, the sequence similarity among KIR genes, the presence of retroviral elements and the dense packing of genes in the KIR cluster might also contribute to the seemingly high recombination activity.²⁶ Other molecular features that could modulate recombination events involve the local GC content and the distribution of DNA-binding motifs for PR domain-containing protein 9 (PRDM9).^28,29 This latter polymorphic zinc finger protein associates with meiotic recombination sites in mammals, contributing to the initiation of double-stranded breaks and thereby facilitating the exchange of genetic material. Allelic variation in PRDM9 concentrates in the DNA-binding domains, resulting in distinct DNA recognition specificities among individuals and species, ultimately shaping differential recombination landscapes.³⁰

In this study, we explored the recombination intervals within the human and macaque KIR regions through pairwise sequence alignment of hybrid genes and recombined haplotypes. By intersecting these recombination hotspots with genomic features that potentially drive rearrangements, we mapped the factors that contribute to the differential diversification observed in the rhesus macaque and human KIR regions.

Materials and methods

DNA extraction, sequencing, and haplotype assembly

We isolated high-molecular-weight and ultra-high-molecular-weight DNA from peripheral blood mononuclear cell samples (±7 × 10⁶ cells) of 13 rhesus macaques of Indian origin that are housed at the Biomedical Primate Research Centre (Rijswijk, the Netherlands) using the Nanobind UL library prep kit (Circulomics). A computational enrichment approach was applied, referred to as adaptive sampling, in concert with sequencing on R9.4 flow cells and an Oxford Nanopore Technologies GridION device. Basecalling was performed with Guppy V3.4.1 software, and the generated reads were filtered on length (>3 kb) and quality score (q-score > 7) before they were mapped to a reference database containing KIR exon sequences from the IPD-KIR database (www.ebi.ac.uk/ipd/kir/) using minimap2. Mapped reads that were generated after adaptive sampling enrichment had an average length 55 kb and were used to construct complete KIR haplotypes by de novo assembly and by alignment to genomic reference sequences using minimap2, reaching an average coverage of 40×. All KIR alleles were annotated based on exon reference sequences or were designated as novel entities based on phylogenetic trees.

Gene density

The gene density of the human and rhesus macaque KIR, LILR, and MHC-DRB regions were calculated for different haplotypes. The regions were defined from the first coding exon on the haplotype to the last coding exon. The length of each region was determined in combination with the number of functional genes and pseudogenes present. For comparison, the gene/length ratios were extrapolated to the number of genes per Mb sequence. Gene densities that were calculated for the different haplotypes were used as independent samples in a t test using RStudio v2024.12.0.³¹ The overall gene densities of the human and rhesus macaque genomes were determined using the total length and annotated genes documented for the reference genomes, GRCH38 and Mmul_10.^32,33 The calculated overall gene density for the rhesus macaque genome might be underestimated due to less comprehensive annotation of the reference sequences compared with the human reference.

Gene similarity

The sequences of KIR genes and their corresponding intergenic stretches were extracted from the KIR haplotypes and aligned using MAFFT.³⁴ Neighbor-joining trees were generated and plotted for the KIR genes and intergenic regions utilizing Geneious Prime software (version 2023.2.1), applying the Jukes-Cantor and the Kimura substitution models. Both models demonstrated similar clustering, substantiated by high bootstrap values. In addition, a python application, SimPlot++, was used to determine the global and local sequence similarity, applying the Jukes-Cantor distance model, a window length of 200 bp, and steps of 20 bp.³⁵ Similarity networks were plotted using the same tool, with the global similarity threshold set to 95% and the local similarity threshold to 100%, analyzing the complete sequence range. The sequence similarity tables were extracted and used to calculate the number of 100% identical intervals between pairs of KIR genes.

DNA-binding predictions for the PRDM9 zinc finger arrays in macaques

Sequence variation of exon 10, encoding the PRDM9 zinc fingers, was investigated for 55 rhesus macaques of Indian origin by amplicon sequencing. Forty of these individuals represented 14 sets of relatives, with at least 1 parent with 1 child, allowing segregation analysis to confirm novel sequences, whereas 15 others were unrelated individuals. Genomic DNA samples of these animals were available from the in-house biobank. Exon 10 was amplified using sequence specific primers and an optimized polymerase chain reaction protocol (Table S1). The amplicons were size selected using gel electrophoresis and extracted from gel using a GeneJet Gel extraction kit (Invitrogen). Tagged amplicons were pooled and purified twice using AMPure XP beads (Beckman Coulter) at a 1:1 bead-to-DNA volume ratio. SMRTbell libraries were generated according to the PacBio Amplicon Template Preparation protocol for circular consensus sequences and sequenced on a PacBio Sequel II platform with P6-C4 sequencing chemistry. An average of 2,000 CCS reads were obtained per individual, which were mapped to a PRDM9 reference sequence extracted from Mmul_10. Novel PRDM9 exon 10 sequences were confirmed by at least 2 independent polymerase chain reactions or by allele segregation in relatives. Eleven sequences encoding the zinc finger array were confirmed, all representing novel PRDM9 alleles (Fig. S1). The new alleles are designated as PRDM9*New, followed by a sequential numbering (e.g. PRDM9*New1). The nucleotide sequences were translated to amino acids for further analysis.

An online Cys₂His₂ zinc finger predictor was used to identify the zinc fingers encoded by the different PRDM9 alleles (HMMER bit scores > 17.7).^36,37 The residues 1, 3, and 6 are reported to facilitate DNA-binding interactions for each zinc finger. The polynomial support vector machine algorithm was applied to predict the DNA-binding motifs for each zinc finger array, for which position weight matrixes were generated. The position weight matrixes were converted to the MEME-motif format, which were subsequently used to generate DNA sequence logos.³⁸ To compare the different motifs, we first aligned them and subsequently stacked the motifs in a phylogenetic tree using the motifStack package in RStudio v2024.12.0 (R Foundation for Statistical Computing).³⁹

Determination of recombination intervals

To determine the site of recombination, sequences of hybrid KIR genes and gene tandems were aligned with their putative donating gene sequences using MAFFT.³⁴ Subsequently, the alignments were used for recombination analysis utilizing RDP4 (Recombination Detection Program version 4) with a sliding window of 30 bp.⁴⁰ The start and end coordinates of the recombination intervals were determined with a 99% confidence interval.

Intersecting repeat elements and PRDM9 motifs with the recombination intervals

We used RepeatMasker (v4.1.4) utilizing the rmBLAST search engine and default settings to screen the generated consensus sequences of each KIR haplotype for repeats.⁴¹ The repeats were categorized into 3 classes: long interspersed elements (LINEs), such as the L1MA4 lineage; short interspersed elements (SINEs), like the AluSx subfamily; and long terminal repeats (LTRs), including MLT1D sequences. In addition, low-complexity stretches and simple repeats were also determined using the same tool. The repeat elements identified by RepeatMasker were listed into a BED file, containing information on the KIR haplotype ID, the coordinates of the repeat, and the type of repeat.

The occurrences of the DNA-binding motifs corresponding to the different PRDM9 allotypes were assessed for all assembled KIR haplotypes by utilizing the FIMO tool.⁴² The default settings were used, and the scan output was converted to a BED file, containing the KIR haplotype ID, the coordinates of the DNA-binding motif, and the ID of the corresponding PRDM9 allotype. The matches with the highest confidence levels (P value <10 × 10⁻⁵ for rhesus macaque motifs and P-value <10 × 10⁻⁶ for human motifs) were split into a second BED file.

The coordinates of the repeat elements and the PRDM9 binding motifs were intersected with the coordinates of the recombination intervals on the KIR haplotypes using BEDTools.⁴³ The “closest” flag was used to identify the features overlapping with the recombination interval, or, when no feature was overlapping, the nearest feature was determined.

Scanning haplotypes for recombination-associated features

Based on the intersections with recombination intervals, 3 LTR elements and 2 broad groups of PRDM9 binding motifs were determined to represent recombination-associated features in rhesus macaques. The 3 LTR elements, 1 within intron 3 and 2 within intron 6, were aligned and plotted in a phylogenetic tree, applying the Jukes-Cantor and Kimura substitution models. Both models demonstrated similar clustering, substantiated by high bootstrap values. The sequence alignment and phylogenetic tree were combined into an illustration using the ggtree package in R.⁴⁴ The 2 PRDM9 motifs groups represent predicted binding specificities of 6 allotypes, including the 2 predominant variants in our characterized rhesus macaque cohort. In addition to these features, sequence motifs were determined that were present in all or in a selection of the recombination intervals, referred to as interval motifs. This motif discovery was performed by utilizing the MEME tool, either forcing the tool to discover motifs shared by all intervals (OOPS flag) or by a selection of the intervals (ZOOPS flag).⁴⁵ The interval motifs with high confidence (P value < 1 × 10⁻²⁰) were selected, aligned, and plotted into a phylogenetic tree as described previously.

The distribution of the recombination-associated features was assessed for all assembled KIR haplotypes. The LTR sequences identified in rhesus macaques were aligned to the haplotypes using minimap2, with the “asm5” flag to only map at a threshold of 95% sequence similarity. The coordinates of the LTRs were extracted from the BAM files and saved as a BED file. In addition, the rhesus macaque and human KIR haplotypes were also scanned for the 2 groups of PRDM9 motifs and the 7 human PRDM9 motifs reported by Altemose et al.,⁴⁶ respectively, using the FIMO tool as described above. The rhesus macaque KIR haplotypes were scanned for the 2 sets of interval motifs as well. To uncover the distribution, the number of occurrences was determined per feature and per KIR gene. This distribution was visualized with a custom script that utilizes the different generated BED files and the circlize package in R.⁴⁷

Results

Characterization of rhesus macaque KIR haplotypes using long read sequencing

The genomic organization of various human KIR haplotypes is well documented,^{14,16,48–50} whereas the characterization of the rhesus macaque equivalent remains relatively limited.²³ Initially, we have deduced over 100 distinct rhesus macaque KIR haplotypes utilizing full-length transcript sequencing in concert with segregation analysis.^18,19,24 This approach revealed the presence of variable numbers and combinations of transcribed KIR genes per deduced region configuration, in addition to allelic polymorphism. We selected 13 of these deduced rhesus macaque haplotypes that display strong indications of recombination events for further analysis. As most macaques are heterozygous for their KIR haplotypes, 3 additional KIR regions with putative standard organizations complemented the selection. The detailed genetic organization of these 16 distinct genomic KIR regions was resolved using long-read sequencing on an ONT platform in combination with adaptive sampling (Fig. 1). In addition, we included 5 KIR haplotypes that we characterized utilizing a Cas9-mediated enrichment protocol.²³

As compared with humans, the rhesus macaque KIR haplotypes display a condensed centromeric and a variably expanded telomeric region, which are separated by a noncoding stretch of DNA. Two framework genes, KIR3DL20 and the pseudogene KIR2DP, are present on nearly all centromeric regions, except for configurations where these loci have been involved in recombination events (H14, H16, H19). The only macaque KIR gene that shares an apparent ortholog with humans is KIR2DL04 and is encountered on most configurations at the start of the telomeric segment. The shortest haplotype contains 5 KIR genes (H4-C) and span 150 kb, whereas the largest region comprises 17 genes (H107) and measures 243 kb. Except for the conserved KIR2DP gene, no other pseudogenes are identified, which indicates a well-orchestrated and highly efficient recombination machinery that diversifies the KIR gene cluster. This contrasts, for example, with the macaque MHC class II region that has many pseudogenes and truncated remnants.⁵¹

To elucidate the molecular mechanisms driving diversification, we set out to identify potential recombination hotspots within the panel of rhesus macaque KIR haplotypes and 14 distinct human KIR region configurations, representing group A and B haplotypes that were previously published by our team and other groups (Fig. 2).^23,48,49

Primate KIR haplotypes display markers of recombination events

Meiotic recombination events that have occurred within the KIR region are characterized by 3 distinct haplotype features: (1) the presence of hybrid KIR genes, in which a donor and acceptor gene were fused in frame; (2) the presence of multiple copies of a specific KIR gene; and (3) the emergence of novel haplotype organizations (Fig. 1). The first haplotype feature marks meiotic rearrangements within genes, whereas the other features may also originate from recombination events in intergenic sequences. In our panel of rhesus macaque KIR haplotypes, we identified 11 hybrid KIR genes, originating from distinct fusion events, distributed over 10 haplotypes (Fig. 1 and Table 1). Three haplotypes, H107, H121, and H130, contain 2 hybrid KIR genes. The generation of these hybrid entities involved 12 different donor genes, including pseudogene Mamu-KIR2DP. Most frequently involved in intragenic recombination events are KIR3DL05 and KIR3DL08, suggesting the presence of strong motifs that promote recombination within their sequences. On 7 haplotypes the presence of multiple copies of 1 or more similar KIR genes indicates the introduction of genes by recombination events. For example, haplotype H9 has 3 copies of KIR3DL07 and 2 copies of KIR3DL05 and KIR3DS02, respectively. We hypothesize that this region configuration resulted from different successive chromosomal rearrangements. Copies of KIR3DL07 and 3DS02 are frequently involved in the expansion of region configurations, again suggesting the presence of strong recombination motifs. Three macaque KIR haplotypes display a peculiar region organization (H14, H16, H19), as their centromeric and telomeric regions are fused, lacking the large noncoding stretch that generally separates KIR2DP and KIR2DL04 (Fig. 1).

The 3 markers of recombination events were also recorded on the selected human KIR haplotypes (Fig. 2). Four distinct hybrid KIR genes were present on 5 region configurations (Table 2), each of them displaying a reduced gene content (Fig. 2). Three of those hybrid KIR entities emerged from reshuffling the segments encoding the binding and signaling domains, while the hybrid KIR2DS1 gene experienced rearrangement only in the exons encoding the leader peptide, originating from KIR2DL1. Several human region configurations containing multiple copies of particular KIR genes have been deduced using standard typing methods in combination with segregation studies,¹⁵ but only 1 of these has been completely characterized at the genomic DNA level (cA01-tB04) (Fig. 2). This expanded region configuration is characterized by the introduction of an entire telomeric segment (tB01). The other selected human KIR region configurations represent standard organizations (groups A and B) or display minor alterations with regard to gene content.

High KIR gene density and sequence similarity may impact recombination rates

In most organisms, an elevated recombination rate is positively correlated with gene density.^52–54 The architecture of the KIR region is characterized by a head-to-tail organization in which densely packed genes are separated only by relatively short intergenic stretches. The studied rhesus macaque KIR haplotypes display a significantly higher gene density than their human counterparts (Fig. 3). This difference is mostly due to longer sequences encountered for introns 5 and 6 in the human KIR genes, which extends the haplotype length. In comparison with other segmentally duplicated immune regions, such as the genes encoding the LILRs (leukocyte immunoglobulin-like receptors), located adjacent to the KIR cluster, and the MHC-DRB genes, the physical proximity of genes is 2-fold higher in the KIR regions of both humans and macaques (Fig. 3). The impact of the significantly denser arrangement of rhesus macaque KIR genes, in contrast to their human equivalents, on potentially increasing the recombination activity remains an enigma.

Sequence similarity and recombination rates are also interconnected factors that modulate the genomic landscape.⁵⁵ The KIR genes and their associated intergenic sequences cluster independently in phylogenetic trees (Fig. S2), reflecting a complex evolutionary history as evidenced by distinctive structural diversity (indels) and point mutations. Beyond these variations, numerous short and long intervals demonstrate a high degree of sequence similarity across KIR genes in rhesus macaques and humans, encompassing exons, introns, and intergenic regions (Fig. S3 and Tables S2 and S3). These highly similar stretches may facilitate unequal alignment of homologous chromosomes during meiosis, thereby enabling the rearrangement of genetic material. In humans, for example, the intergenic sequences associated with KIR2DL5B and 3DP1 are highly similar (Fig. S4), which may facilitate recombination events that reduce the gene content of the centromeric cB02 region (Fig. 2). Similarly, in rhesus macaques, KIR3DS02 is often involved in rearrangements and shares at least 1 identical stretch of 200 bp with 13 different KIR genes. In contrast, KIR3DL20 and KIR2DL04 are stable framework genes on most macaque haplotypes and share identical sequences with only 1 other KIR gene, KIR3DS06 and KIR3DS01, respectively. However, these latter genes are not involved in the formation of the hybrid entities that are documented for the framework genes. Additionally, more identical stretches are shared among sets of human KIR genes than among rhesus macaques genes (Table S3), but there are no indications that they often trigger recombination events in the contemporary time. Therefore, it remains uncertain whether the high sequence similarity between KIR genes contributes to local recombination activity.

The architectural features involving gene density and sequence similarity make the KIR regions in primates susceptible to meiotic recombination. However, these features alone do not account for the precisely orchestrated process, which only generates in-frame hybrid KIR genes and intact expanded and contracted haplotype organizations. Instead, it is likely that motifs at or near recombination sites play a more dominant role in guiding the rearrangement of KIR region configurations.

Determination of the recombination sites

Recombination sites, defined as short sequence intervals in which the DNA strand breaks during meiosis, can be identified within the KIR region through the pairwise alignment of rearranged segments and their potential precursor sequences. We mapped the recombination interval for the 11 rhesus macaque hybrid KIR genes in our panel (Table 1). For example, a recombination interval of 72 bp is located in intron 3 of KIR3DLW43, a hybrid entity that consists of segments originating from KIR3DL02 and KIR3DL08 (Fig. 4A). Other intervals range from approximately 35 to 450 bp in length and occur both in introns and exons. The length of recombination intervals increases with the level of sequence similarity between the 2 donating genes, as discriminative SNPs are required to pinpoint putative meiotic breakpoints. Recombination intervals were also determined for the 4 human hybrid KIR genes present in our panel (Table 2).

In addition to the generation of hybrid genes, recombination events are also responsible for the expansion and contraction of KIR region configurations by the introduction and deletion of genes. The breakpoints of these rearrangements are likely situated near or within intergenic sequences. To identify such recombination intervals, we first assessed the occurrence of pairs of neighboring KIR genes. Notably, most KIR genes exhibit strong linkage with their adjacent gene (Table S4). For example, KIR3DL08 is adjacent to KIR3DL01 in all of the haplotypes containing this gene in our macaque panel. Uncommon KIR gene tandems may mark meiotic recombination events, 10 of such were present in our cohort (Table 3). For instance, KIR3DL07 is often linked to the 3′ side of KIR3DL05, but in one case this position was occupied by KIR3DS02 (H69), which may indicate a recombination event. The associated breakpoint probably locates in a recombination interval of 382 bp that was determined for the KIR3DL05–KIR3DS02 tandem in intron 6 of the former gene (Fig. 4B). Most recombination intervals in the other potentially rearranged gene tandems also start within intron 6 of the former gene and may extend into introns 7 and 8. The location of this putative hotspot is supported by pairwise sequence alignments of 3 distinct gene tandems involving KIR2DL04 (Fig. S5). The recombination events mapping to intron 6 may have generated previously unrecognized hybrid KIR genes, in which the exons encoding the signaling domains are swapped. This intragenic process is further substantiated in humans by the tight linkage of intergenic sequences and their associated genes (Fig. S4), suggesting a restriction on the shuffling of KIR promoters. Two of the human hybrid KIR genes, KIR2DL1 (KU645198) and KIR2DS2 (KP420437), also share recombination intervals within intron 6, but this hotspot seems to be less active.

Although precise meiotic breakpoints could not be determined, the identified recombination intervals indicate that all rearrangements in the macaque KIR region generate hybrid gene entities. Aligning sequence features that may promote recombination with the defined intervals might pinpoint the elements that contribute to the dynamic nature of the macaque and human KIR regions.

DNA-binding motifs of rhesus macaque and human PRDM9 allotypes

An important modulator of meiotic recombination events, PRDM9, exhibits allelic variation in its zinc finger array, which may affect its DNA recognition motif.⁵⁶ So far, only a single rhesus macaque PRDM9 sequence is published, and the extent of allelic variability remains to be explored. We investigated the genetic diversity of the PRDM9 zinc finger array in 55 Indian rhesus macaques and identified 11 novel alleles, 2 of which were predominantly present in our cohort (Fig. 5). The nucleotide variation is mostly confined to the stretches encoding the different zinc fingers, and in particular concentrates around the 3 DNA-contacting residues (Fig. S1). The discovered alleles encode for 8 distinct zinc finger arrays, and their corresponding DNA-binding motifs were predicted computationally (Fig. 6). This revealed 4 broad motif groups predicted to be binding targets for the different PRDM9 allotypes detected in rhesus macaques. The 2 most frequent alleles in our cohort, New1 and New2, have distinct predicted binding specificities.

In humans, over 40 PRDM9 alleles are recorded.^57,58 The DNA recognition motifs of the 3 most prevalent alleles, designated as PRDM9-A, PRDM9-B, and PRDM9-C, are extensively examined.^46,57 The A and B alleles broadly recognize the same 13-mer sequence and are predominant in non-African populations. Within African cohorts, the PRDM9-A allele manifests at a frequency of 50%, complemented by 12% prevalence of PRDM9-C and a wide diversity of alleles with lower frequencies.⁵⁹ Seven motifs with a close internal match to the 13-mer sequence are recognized by the PRDM9-B allele.⁴⁶ These binding motifs are likely shared with PRDM9-A, which only differs by 1 DNA-contacting residue in its zinc finger array.

The PRDM9 binding domains in rhesus macaques may contain 9 to 11 zinc fingers, whereas the 3 prevalent human alleles exhibit 13 or 14 zinc fingers (Fig. 5). In addition, except for the proximal zinc finger that is not involved in DNA-binding, only 1 zinc finger (DSY) is conserved in the arrays of rhesus macaques and humans, sharing identical DNA-contacting residues. These observations support the distinct binding motifs that are determined (or predicted) for their PRDM9 alleles (Fig. 6). Consequently, different recombination landscapes are expected for rhesus macaques and humans, which is in line with a comparative study highlighting a variable distribution of recombination hotspots.⁶⁰ This suggests that allelic variability of PRDM9 may impact the differential diversification of KIR haplotypes at a species level.

Repeat elements and PRDM9 motifs associate with recombination intervals

A total of 21 recombination intervals were identified in the panel of rhesus macaque KIR haplotypes, whereas 4 intervals were determined in the selected human equivalents. For each interval, we assessed the proximity to, or overlap with, repeat elements and PRDM9 motifs.

Twelve macaque intervals intersected with repeat elements, including all recombination sites within intron 6, except for 1 (Tables 1 and 3). Predominantly, the repeat elements represent LTRs, complemented by 2 LINE and 1 SINE repeats. Two distinct LTR sequences starting in intron 6 comprise a relatively short stretch of approximately 240 bp and a longer one extending up to 910 bp, respectively (Fig. S6). Both LTRs in intron 6 are different in comparison with the intersecting LTRs present in intron 3. An example of this latter LTR maps within the interval determined for the hybrid KIR3DLW43 gene (Fig. 7). If not overlapping, the nearest repeat elements are situated at an average distance of 551 bp from the recombination intervals, with proximities ranging from 64 to 1,182 bp. Notably, all of these proximal elements belong to the LTR family, with a single exception that belongs to the SINE family (Table 1). Three of the 4 human recombination intervals intersected with a repeat element from the SINE, LINE, or LTR class (Table 2). The largest interval aligned with 2 repeat elements, a SINE and an LTR. The LTR found in intron 6 of the hybrid KIR2DS2*005 gene exhibits a sequence similarity of up to 88% with its counterparts in rhesus macaques, indicating a retroviral element that integrated long before speciation. These type of ancient LTR integrations are more often associated with the diversification of immune regions, as has been demonstrated for the MHC-DRB genes in rhesus macaques.⁵¹

Eight recombination intervals in rhesus macaques intersect with a predicted PRDM9 motif, 5 of which also align with a repeat element (Tables 1 and 3). Most intersecting binding motifs correspond to PRDM9*New7, PRDM9*New10, and PRDM9*New9 (Fig. 6). When no motif overlapped with a recombination interval, we identified the nearest motif at an average distance of 200 bp, with proximities ranging from 22 to 542 bp. For example, a motif corresponding to PRDM9*New2 was identified 95 bp upstream of the recombination interval of KIR3DLW43 (Table 1 and Fig. 7). Hence, PRDM9 typically does not bind directly at the exact site where double-stranded breaks form during meiosis. Instead, it commonly binds to DNA sequences nearby, usually within a few hundred base pairs from the double-stranded break site, as was previously documented in humans.^60,61 Consequently, also the proximal motifs that do not directly overlap with a recombination interval might contribute to a recombination event. Sets of similar PRDM9 binding motifs, categorized into motif groups 1 and 2, were primarily found near or overlapping with the 21 recombination intervals (Tables 1 and 3 and Fig. 6). In contrast, the motifs associated with the New3 and New4 allotypes were only identified once near an interval. The predominant motif groups both associate with 1 of the 2 most frequently identified PRDM9 alleles (New1 and New2) present in our cohort of macaques. The DNA-binding motifs of human PRDM9 intersected with 2 recombination intervals, whereas for the other 2 putative rearrangements a motif was identified at close proximity (Table 2). Given the slight differences among the human PRDM9 binding specificities, multiple motifs intersected with 2 of the intervals.

When neither a repeat element nor a PRDM9 recognition motif overlapped with a recombination interval in rhesus macaques, which was the case for 6 intervals, the nearest feature was identified at a minimal distance ranging from 22 to 183 bp. This mostly involved a proximal PRDM9 binding motif.

Predicting putative recombination sites on KIR haplotypes

The characterization of complete KIR haplotypes facilitated the determination of recombination intervals and the identification of diverse molecular features associated with them. By mapping these features on the genomically characterized KIR haplotypes, we could assess the presence of other potential breakpoints that may initiate recombination during meiosis. In addition, we determined 2 sequence motifs that recurred in all 21 recombination intervals in rhesus macaques, as well as 2 motifs that were recurring in at least 6 of the intervals (Fig. S7). The low conservation across all recombination intervals suggests that diversification of the KIR region is not driven by a single molecular feature. Nevertheless, these recurring motifs might serve as indicators for the presence of potential recombination hotspots on the KIR haplotypes.

Mapping LTR sequences, PRDM9 binding motifs, and recurring interval motifs onto the KIR haplotypes revealed variable distributions of recombination-associated features at the gene level, exemplified by rhesus macaque KIR haplotype H5 (Fig. 8 and Table 4). Overall, the LTR element intersecting with recombination intervals in intron 3 was found in 9 distinct KIR genes (Table 4). Other genes may have an LTR element at this position, but these are distinct from the sequence that was identified to intersect with recombination intervals (<95% sequence similarity). In our macaque panel, 3 hybrid KIR genes originate from a rearrangement with a double-stranded break in intron 3 (Table 1), and their generation may have been facilitated by the presence of these LTR elements. Other recombination-promoting LTR elements are situated in intron 6 and are present in the majority of KIR genes, appearing as either short or long variants. These highly similar repeats possibly contribute to the continuous expansion and contraction of KIR haplotypes by promoting the reshuffling of genes through recombination. This observation is substantiated by the highly conserved organization of the centromeric haplotype, containing KIR3DL20 and KIR1D, which both lack the recombination-associated LTR element in intron 6. Furthermore, KIR3DL08, also lacking the LTR in intron 6, consistently appears adjacent to the 5′ site of KIR3DL01 (Table S4). This suggests that any gene lacking this repeat element may not participate in rearrangements. Approximately two-thirds of the macaque KIR genes possesses at least 1 PRDM9 motif with high confidence, mostly involving recognition patterns of the PRDM9*New2, PRDM9*New6, and PRDM9*New9 (group 2) allotypes. Most recognition motifs are located in introns 4 and 6 and thereby may contribute to the formation of hybrid KIR genes and the expansion and contraction of haplotypes. This is supported by the observation that most PRDM9 motifs are located in KIR3DS02, which is indeed frequently involved in the expansion of haplotypes. Predicting the potential involvement of PRDM9 motifs identified with lower confidence in initiating recombination events is difficult. However, given that some of these motifs were detected near recombination intervals makes it is plausible that they may play an active role as well. The motifs that were generated from the 21 recombination intervals in rhesus macaques, or a selection thereof, were identified in most KIR genes. These motifs mostly mapped to exon 4, the end of intron 5, and within intron 6, predicting potential recombination hotspots that are supported by several detected recombination events.

Several sets of KIR genes display patterns of co-occurrence and co-exclusion.²⁴ A strong linkage is, for example, documented for KIR3DS05, KIR3DL10, and KIR3DL01 in macaques, indicating the absence of recombination events within this cluster (Fig. 1). A recombination-associated LTR was not identified in intron 3 of these linked genes. However, in intron 6, either a short or long LTR was detected. Furthermore, while PRDM9 motifs are absent in KIR3DL10, 1 and 2 of these motifs were identified in KIR3DS05 and KIR3DL01, respectively. Other genes are involved in similarly linked clusters, like KIR3DS02, KIR3DL05, and KIR3DL07, with again the central gene lacking PRDM9 motifs, while different LTR sequences are found present in all genes. These observations suggests that not all recombination-associated elements may be active, or that selection avoids the uncoupling of these genes.

Equivalents of the LTR in intron 3 were only identified in human KIR3DL2, whereas the LTR elements associated with intron 6 were found in 7 distinct human KIR genes (Table 5). However, the rhesus macaque and human LTR elements deviate in their sequence length, which might have an impact on the initiation of recombination.⁶² A PRDM9 motif was identified with high confidence in all human KIR genes, except for KIR2DL5 and KIR2DS4. However, due to the limited number of recombination intervals determined in the available human KIR haplotypes, establishing a correlation with these molecular features is at present challenging.

Overall, it appears that recombination in the KIR region is driven by a combination of molecular features, including specific LTR elements and PRDM9 binding motifs that may play a dominant role.

Discussion

The thorough characterization of complex immune regions, such as the KIR region, was always a challenging enterprise using short-read sequencing platforms, due to copy number variation, high levels of sequence similarity, and extensive allelic polymorphism. The most common technique to characterize the KIR repertoire in humans is based on the determination of the presence or absence of gene segments using sequence-specific primers, which has its limitations.⁶³ More in-depth studies in other primate species involved partial or full-length KIR transcript sequencing, occasionally in concert with segregation studies to resolve haplotype information.¹⁸ This approach revealed different levels of complexity, but precise KIR haplotype organizations at the genomic level remained an enigma. Long-read sequencing strategies opened avenues to assess and phase complex immune regions. We applied, for instance, a Cas9-mediated enrichment approach in combination with ONT sequencing to assemble phased rhesus macaque and human KIR haplotypes.²³ Although considerable progress was made, this strategy is relatively expensive and time-consuming. In this communication a rapid software-controlled enrichment strategy is presented, known as adaptive sampling, that successfully resolved completely phased KIR haplotypes in a panel of rhesus macaques. These animals have an expanded MHC class I region as compared with humans.²⁵ We aimed to determine if this expansion impacts the extent and evolution of the KIR repertoire. On average, the evolution of the human KIR region seems to be under less selective pressure as compared with rhesus macaques. Moreover, the growing number of completely characterized KIR haplotypes offers the opportunity to explore the molecular features that drive the expansion and contraction of this immune region.

The elevated recombination activity observed in the KIR regions of rhesus macaques and humans, compared with other segmentally duplicated immune gene systems, may involve their high levels of gene density and sequence similarity (Fig. 3; Fig. S3). However, it is more likely that additionally an even more sophisticated mechanism is at play, which maintains functional stability by producing in-frame recombination events, while generating genetic diversity. This mechanism appears to involve the presence of PRDM9-binding sites and LTR elements, which are located within introns, while the exon boundaries remain conserved. These features may shape the recombination landscape of KIR genes and haplotypes by facilitating the exchange of genomic material at specific hotspots during meiosis.

However, not all of these recombination-associated elements are consistently active. For example, PRDM9 is expressed in germ cells during meiosis, but only a fraction of the recombination hotspots are activated per cell. Moreover, the binding of PRDM9 alone does not guarantee the occurrence of a double-stranded break, as the interplay of various factors within the recombination machinery is essential for facilitating successful rearrangement.⁶⁴ Active binding sites of PRDM9 in gametes may be identified by a parallel detection of DMC1 binding, a meiosis-specific recombinase that catalyzes strand invasion necessary for recombination.⁶⁵ While this combined information is currently unavailable for the KIR regions in humans or macaques, future studies employing this approach could offer valuable insights into the active PRDM9 hotspots that drive diversification of KIR haplotypes. In addition, this exercise would also provide a more accurate pattern of the PRDM9 binding motifs. In this communication, we applied a widely-used and validated prediction tool to define the binding motifs of the various PRDM9 alleles identified in the panel of rhesus macaques.³⁶ This software tool correctly predicted the binding of the most prevalent human variant, PRDM9-A, to a 13-mer motif that is enriched at recombination hotspots.³⁰ Although these in silico predictions offer indicative binding motifs, their reliability may be limited due to the complexity of the interaction between the diverse array of zinc fingers and DNA. This is exemplified by the diverse set of DNA-binding motifs identified for human PRDM9-B, sharing a close match to the 13-mer motif but displaying distinct internal spacings.⁴⁶ The diverse binding specificities may be explained by different factors, including the varied combinations of zinc fingers and the potential formation of PRDM9 multimers.^66,67 The binding specificities of rhesus macaque and human PRDM9 are clearly distinct, stemming not only from differences in the type of zinc fingers but also in the number of fingers within each array. The smaller array in rhesus macaques correlates with a narrower DNA-binding motif, potentially leading to its increased occurrence compared with the longer motifs associated with human PRDM9. Overall, the overlap of PRDM9 binding targets with, or in proximity to, all 21 examined recombination intervals in the KIR region of rhesus macaques supports the role of this zinc finger protein in facilitating recombination events.

The other feature implicated in modulating the diversification of the KIR region is the distribution of specific LTR elements. These retroviral elements appear to be predominantly involved in recombination events within introns 3 and 6 of most KIR genes, represented by either short or long LTR sequences. Recently, we reported a similar LTR-driven recombination mechanism in the MHC-DRB region of rhesus macaques.⁵¹ In this context, 2 relatively long LTR sequences within and adjacent to the pseudogene DRB6 were associated with the shuffling of sets of paralogous DRB genes. In contrast to the KIR region, a relatively high number of inactive gene copies are generated by the diversification of the MHC-DRB cluster. This discrepancy indicates that a more sophisticated process diversifies the KIR region. Variation in the lengths and distributions of LTR elements in these 2 immune regions might contribute to the differential recombination events.

In addition to variation in molecular features that may promote rearrangements, the continuous diversification of the KIR gene content in macaques, compared with the more confined haplotypes in humans, is likely modulated by evolutionary selection. Selective pressure favoring functional KIR genes and haplotypes would probably drive the elimination of nonproductive rearrangements from the population by natural selection. Consequently, out-of-frame recombination events may occur but are less likely to be detected. However, an example of a dysfunctional gene that originated from a rearrangement and remained in human populations is KIR2DL5B*002.^68,69 Although this gene is in-frame, it is not expressed due to the introduction of the promoter region from pseudogene KIR3DP1. This suggests that inactive alleles can occasionally arise from recombination events and remain in the gene repertoire, even though this occurrence is rare.

Different selective pressures may have an effect on the KIR gene repertoire. In humans, the selection pressures may largely involve successful pregnancy, as the KIRs play a pivotal role in the maternal immune tolerance.^70,71 Specifically, KIRs bind to certain epitopes on HLA-C molecules, which are preferentially expressed on trophoblasts. This interaction facilitates processes like trophoblast invasion and vascular remodeling, both essential for a successful pregnancy. The evolutionary pressure to optimize these interactions has led to the selection of specific KIR haplotypes that promote reproductive success, resulting eventually in a less diverse but more specialized KIR repertoire in humans. In contrast, the KIR system in rhesus macaques is primarily driven by the need to maintain a functional relationship with the expanded MHC-A and MHC-B molecules, thereby indirectly adapting to the rapid evolution of pathogens, including an array of retroviruses.^5,72 This difference in primary selective pressures is further enlarged by the education pathways of NK cells. Unlike humans, where KIR education is heavily reliant on HLA-C interactions, rhesus macaque NK cells are primarily educated through the NKG2A pathway interacting with MHC-E molecules.^72,73 The shift in NK cell education pathways may reduce the selective pressure on specific KIR-MHC interactions in macaques. This might allow for greater expansion and diversification of the KIR gene repertoire without the reproductive constraints seen in humans.

Conclusively, the expansion and contraction of the primate KIR region appears to enhance the evolutionary adaptability, while the underlying molecular process and selective pressure maintain the functional integrity of KIR genes. We identified sequence features that may collectively contribute to the sophisticated diversification within the KIR region of primate species, including specific LTR elements and sets of PRDM9 binding motifs. Variation in these sequences across species may account for the contrasting levels of diversity recorded for the rhesus macaque and human KIR haplotypes.

Acknowledgements

The authors thank F. van Hassel for design of figures and artwork.

Author contributions

Marit K.H. van der Wiel (Formal analysis [equal], Investigation [equal], Methodology [equal], Validation [equal]), Ngoc Giang Le (Formal analysis [equal], Software [lead], Visualization [lead]), Nanine de Groot (Investigation [supporting]), Natasja G. de Groot (Conceptualization [supporting]), Ronald E. Bontrop (Conceptualization [supporting]), and Jesse Bruijnesteijn (Conceptualization [lead], Formal analysis [lead], Investigation [equal], Methodology [lead], Project administration [equal], Supervision [lead], Visualization [equal])

Supplementary material

Supplementary material is available at The Journal of Immunology online.

Funding

None declared.

Conflicts of interest

None declared.

Data availability

All sequencing data generated in this study has been submitted to the European Nucleotide Archive (https://www.ebi.ac.uk/ena/browser/home) under project number PRJEB77196. The annotated consensus sequences of the rhesus macaque KIR haplotypes have been submitted under accession numbers ERZ24798745 to ERZ24798765.

References