Human and Rhesus Macaque KIR Haplotypes Defined by Their Transcriptomes

Killer-cell Ig-like receptors (KIRs) are expressed on NK cells and subsets of T cells (1, 2), and play a key role in immune recognition through interactions with the highly polymorphic MHC class I molecules (3, 4). For example, KIRs may play an important role in the detection of aberrant MHC class I expression on tumor and virally infected cells, and their subsequent elimination (5, 6). KIRs are type I transmembrane glycoproteins that consist of two or three extracellular Ig-like domains (2D or 3D) as well as a stem, transmembrane region, and cytoplasmic tail. The length of the cytoplasmic tail can be either long (L), including two ITIM motifs, or short (S), and characterizes inhibitory or activating KIRs, respectively. KIRs with one extracellular Ig-like domain and a truncated cytoplasmic tail (KIR1D) are observed in some nonhuman primate species, and seem to have no counterpart in humans (7, 8). The nomenclature of the KIR genes is based on the functional (L or S) and structural (2D or 3D) characteristics, and takes into account allelic variation as well (9).

In humans, the KIR gene cluster is located within the leukocyte receptor complex on chromosome 19q13.4, and displays copy number variation (CNV) at the population level, as reflected by a variable number of tandemly arranged KIR genes (10, 11). A KIR haplotype contains 7–12 expressed genes, three of which are considered framework genes: KIR3DL3, KIR3DL2, and KIR2DL4 (12). A fourth framework gene is KIR3DP1, which is a pseudogene. Based on the genetic makeup, the human KIR haplotypes can be categorized into two groups (13). Group A haplotypes are characterized by seven KIR genes, including the framework genes and only the activating KIR2DS4 structure, whereas group B haplotypes can contain up to 12 genes including the framework genes and multiple activating receptors. The KIR genes can be further divided into phylogenetic lineages (I, II, III, and V), each characterized by structure and MHC class I specificity; lineage I includes KIR2DL4 and KIR2DL5, lineage II includes KIR3DL1/S1 and KIR3DL2, lineage III includes KIR2DS1/2, KIR2DL1/2/3, KIR2DS3/5, KIR2DS4, and the pseudogenes, and lineage V includes KIR3DL3.

In addition to the CNV, allelic polymorphism is another important feature of the KIR gene system. In humans, the greatest expansion is observed for KIR lineage III genes in both the telomeric and centromeric region of the haplotype. A total of 15 human KIR genes and two pseudogenes have been characterized, and up to 907 unique full-length KIR alleles have been cataloged (14).

As various KIRs may bind specific but differential structures on MHC class I molecules, the KIR repertoire influences, in part, the variability of the immune response. Because both the MHC and KIR gene systems display substantial levels of polymorphism and segregate as independent entities located on different chromosomes, the potential repertoire of MHC-KIR interactions may vary considerably, even among related individuals within a family. The presence or absence of certain KIR alleles and MHC-KIR interactions are associated with disease susceptibility and its progress, but may also play a role in transplantation and reproductive biology (5, 6, 15, 16).

In recent years, our understanding of the biology and evolution of the KIR gene system has greatly expanded, although some key questions remain to be answered. Suitable animal models to study KIR-related diseases are more or less confined to nonhuman primate species, because rodent species have another system executing similar tasks that arose as a result of convergent evolution (17, 18). Macaques, for example, share a close evolutionary relationship with humans, which is evidenced by similar pathology and immune responses in models for infectious and autoimmune diseases (19–21). Initial genomic characterization of the KIR gene repertoire in rhesus macaques (Macaca mulatta, Mamu) highlighted substantial similarities along with some differences as compared with humans (8, 22–24). For example, the macaque KIR gene cluster shows an extreme expansion of lineage II KIR genes, mainly in the telomeric part of the haplotype, which might be associated with the multiplicity of the KIR-interacting MHC-A and -B genes (25, 26). This extensive gene CNV exceeds the lineage III KIR expansion observed in humans. Thus far, 22 KIR genes and 218 alleles have been reported in macaques (7, 22, 27–35). KIR haplotypes in macaques can contain up to 11 genes, some including Mamu-KIR2DL04 and Mamu-KIR1D, which belong to KIR lineages I and III, respectively. The human KIR haplotypic division differentiating between the more activating (B) and the more inhibitory (A) haplotypes is not as obvious for macaques (22, 26).

A limited number of macaque haplotypes has been characterized by studying segregation in families in combination with conventional sequencing methods, leading to haplotype definitions that were based on the presence of both partial and full-length cDNA sequences. Although these methods provided insights, they were not always sufficient to either resolve allele-level haplotypes or to identify genes with low transcription levels. In contrast, most human KIR haplotypes were characterized by determining the presence or absence of KIR genes at the genomic DNA (gDNA) level. As a consequence, crucial information on the allele level, CNV, transcription level, and transcriptional modifications, such as splicing, may be missed. As particular KIR alleles are expected to be associated with health and disease, a comprehensive method is required to characterize the complete KIR transcriptome. In this study, we report a next-generation single-molecule real-time (SMRT) sequencing method on the Pacific Biosciences (PacBio) Sequel platform, which allowed us to obtain full-length KIR transcriptomes, as well as KIR haplotypes for both human and rhesus macaque families. This approach provides a significant step forward, which may aid in a better understanding and interpretation of KIR allele–associated diseases, as well as the immune reactivity in transplantation and reproductive biology.

A large pedigree-based Indian rhesus macaque family, with a total of 30 animals, was selected from the self-sustaining colony housed at the Biomedical Primate Research Centre (Fig. 1). EDTA or heparin whole blood samples were obtained during regular annual health checks, and PBMCs were isolated from heparin blood samples. PBMCs of 15 related humans were provided by the immunohematology and blood transfusion department of the Leiden University Medical Center (Fig. 4). Informed consent was obtained from all participants.

Total RNA was extracted directly from EDTA whole blood samples or from ± 15 × 10⁶ PBMCs with RNeasy Mini Kit (Qiagen, Valencia, CA) in accordance with the manufacturer’s instructions. First-strand cDNA was synthesized with the RevertAid First Strand cDNA Synthesis Kit (Invitrogen, Carlsbad, CA) using oligo(dT)₁₈ primers. gDNA was extracted from EDTA whole blood samples by a standard salting-out procedure, or from ± 15 × 10⁶ PBMCs with an AllPrep RNA/DNA Mini Kit (Qiagen) according to the manufacturer’s instructions.

PCR with different primer sets (Table I) was performed on cDNA using Phusion Hot Start II DNA Polymerase (Thermo Fisher Scientific, Waltham, MA) to obtain full-length macaque and human KIR amplicons. Each primer was tagged at the 5′ end with a 16 bp barcode, designed for the PacBio platform, to identify samples by unique barcode combinations. Thermal cycling conditions were denaturation at 98°C for 2 min, followed by 32 cycles of 98°C for 20 s, 66°C for 45 s, and 72°C for 2 min, except for the primer set derived from Moreland et al. (24): denaturation at 98°C for 2 min; five cycles of 98°C for 20 s, 68°C for 5 s, 66°C for 5 s, 63°C for 30 s, 60°C for 5 s, and 72°C for 2 min; 29 cycles of 98°C for 20 s, 63°C for 30 s, and 72°C for 2 min. Appropriately sized PCR products of ∼1250 bp were selected by gel electrophoresis and purified using a GeneJet Gel extraction kit (Invitrogen). The amplified KIR amplicons were pooled and purified twice, using AMPure XP beads (Beckman Coulter, Woerden, the Netherlands) at a 1:1 bead to DNA volume ratio. The DNA concentration of purified pooled samples (>1 μg total DNA) was measured using the Qubit dsDNA HS assay kit and Qubit 2.0 Fluorometer (Thermo Fisher Scientific).

SMRTbell libraries were generated according to the PacBio Amplicon Template Preparation protocol for circular consensus sequences (36). In brief, PCR product pools were end-repaired, and hairpin adapters were subsequently incorporated using the PacBio DNA Template Prep Kit 2.0. After the removal of failed ligation products, the SMRTbell templates were purified with 0.6× volume of AMPure XP beads. The volume of sequencing primer and polymerase was determined using the PacBio calculator. Polymerase-bound SMRTbell libraries were MagBead loaded over zero-mode waveguides and immobilized. SMRTbell library generation and sequencing were performed by the Leiden Genome Technology Center using a PacBio Sequel instrument with P6-C4 sequencing chemistry.

Circular consensus sequences were obtained and selected for high read quality (value of 0.99 or higher). The data were demultiplexed based on unique barcoding, and were used to type individual samples.

Geneious Pro R10 software (37) was used to map the reads to a reference database, including all reported full-length and partial human or macaque KIR cDNA sequences (7, 14, 22, 27–35), to identify 100% matched reads (100% overlap, 0% mismatch, maximum ambiguity = 1). The remaining reads were de novo assembled, and the consensus of each contig was phylogenetically aligned with the human or macaque KIR sequence reference database. Official designations of novel sequences were determined by comparison with the reference sequence databases and by phylogenetic analysis, and were provided by the Immuno Polymorphism Database, which will shortly release a database with nonhuman primate KIR sequences. The novel sequences were confirmed when observed in two independent PCRs or shown to segregate into families, and were subsequently submitted to the European Nucleotide Archive and assigned an accession number (https://www.ebi.ac.uk/ena/).

Alleles that could not be confirmed by segregation and had low PacBio read counts (<3 reads) or alleles that seemed to be duplicated on a single haplotype were confirmed by Sanger sequencing. To amplify products that distinguished the different genes or alleles, specific primers were designed within the exons (Supplemental Table I). For all primers, the PCR conditions were denaturation at 98°C for 2 min, followed by 32 cycles of 98°C for 20 s, 66°C for 45 s, and 72°C for 2 min. PCR products were subjected to gel electrophoresis, and bands of the appropriate size were purified with a GeneJet Gel extraction kit. For Mamu-KIR3DL07, Mamu-KIR3DL05, and Mamu-KIR3DS02 cloning was performed as previously described (28). Sequencing of the PCR products or isolated cloned amplicons was performed on a 3500XL Genetic Analyzer automatic sequencer (Applied Biosystems, Foster City, CA). Sequences were analyzed with SeqMan Pro (DNASTAR, Madison, WI) and MacVector (MacVector, Cambridge, U.K.) software.

In addition to typing at the transcription level, both human and macaque samples were also typed at the gDNA level for the presence or absence of KIR genes. Rhesus macaque DNA samples were typed using quantitative PCR and melt curve analysis as previously described (23). Human KIR genotyping was performed by Sanquin (Department of Immunogenetics, the Netherlands) using the Olerup SSP KIR typing kit (Olerup SSP AB, Stockholm, Sweden) in accordance with the manufacturer’s instructions.

The rhesus macaque KIR system was previously characterized by cloning and Sanger sequencing, Roche/454 pyrosequencing, and microsatellite analysis, which resulted in a database of 218 partial and full-length Mamu-KIR alleles (22–24, 27, 33, 38). Validation of these sequences was provided by independent amplifications and, whenever possible, segregation studies. However, these characterization strategies were time-consuming, and were often insufficient to resolve full-length allele sequences.

Experience taught that the PacBio Sequel platform offers a substantial number of improvements compared with conventional sequencing strategies, such as higher throughput and longer reads with high accuracy by circular consensus sequencing (39, 40). Taking these advantages into account, we set up a pipeline to characterize the KIR gene system using a PacBio Sequel platform. Initially, we calibrated the PacBio platform by reanalyzing macaque blood samples that had been previously typed for KIR by a conventional methodology. With the current PacBio sequencing protocol, we confirmed in considerably less time the rhesus macaque KIR results that had been obtained earlier. Furthermore, additional KIR alleles and genes were identified, and partial sequences that had been missed by conventional methods could be extended. In the following set-ups, we used a family segregation concept, so that identical sequences could be obtained and confirmed by analyzing different but related individuals.

The previously published macaque KIR3D/1D and KIR2DL04 primers located in conserved regions of the 5′ and 3′ untranslated regions (UTRs) were barcoded with PacBio sequence tags, and used to amplify full-length KIR cDNA transcripts (22, 24). A combination of two generic Mamu-KIR3D/1D and one specific Mamu-KIR2DL04 primer sets was required to amplify the complete KIR transcriptome in macaques (Table I).

A family of 30 macaques was selected that originated from eight founders and covered four generations, which allowed an extensive segregation study (Fig. 1). An average of ∼9000 Mamu-KIR3D/1D and 5800 Mamu-KIR2DL04 reads were obtained per animal (±100,000 reads per PacBio Sequel cell), of which 20–45% were perfectly matched with the Mamu-KIR allele library that consisted of 218 annotated Mamu-KIR3D/1D and Mamu-KIR2DL04 sequences. The remaining reads that did not match with the Mamu-KIR library were novel alleles, partial KIR sequences, or KIR sequences containing random single nucleotide gaps, which had been introduced by PacBio sequencing. A total of 29 unreported Mamu-KIR alleles were identified in 30 related rhesus macaques (Tables II, III). Six new KIR alleles were identified for Mamu-KIR3DL07 and Mamu-KIR3DL20, whereas for the other KIR genes, two, one, or no new alleles were detected.

In two of the newly detected alleles, premature stop codons were identified, but both were located in exon 9. For example, an insertion of a single cytosine was found at base pair position 1014 in Mamu-KIR3DL11*009 (LT906600), which resulted in a frame shift and the introduction of an early stop codon at base pair position 1171 (exon 9). This cytosine insertion was confirmed using multiple independent PCRs, and the frame shift might suggest a Mamu-KIR3DL11 isoform with a truncated cytoplasmic tail, thereby lacking the ability to signal via ITIMs. In Mamu-KIR3DL05*013 (LT906588), a point mutation at base pair position 1294 (exon 9) introduced a stop codon, which resulted in a transcript likely to encode a KIR3DL protein with only one ITIM.

In addition to novel alleles, six previously reported Mamu-KIR sequences were extended, and another 30 alleles were confirmed. Overall, of the 69 KIR genes/alleles identified in the rhesus macaque family (Table II), almost half of them were novel, illustrating the power of the platform and the extensive allelic polymorphism of the macaque KIR gene system.

In addition to allelic polymorphism and CNV, the plasticity of the KIR gene system is also reflected by recombination events, such as inter- and intrachromosomal rearrangements, which might result in the formation of recombinant in-frame KIR genes (Fig. 2A, 2B) (41–44). This type of generation of novel gene entities, caused by the fusion of different genes, has been described in humans, but thus far has not been observed in rhesus macaques. In this study, comprising only one extended macaque family, four novel recombinant Mamu-KIR sequences were identified at the transcription level, composed of segments from two different KIR genes (Fig. 2C). Two recombinant genes consisted of the 5′ end to exon 4 of Mamu-KIR3DL10 and exon 5 up to the 3′ end of Mamu-KIR3DL02, but seemed to originate from two independent fusion events. In these two recombinant genes, most point mutations were observed in the Mamu-KIR3DL10 segment, whereas the Mamu-KIR3DL02 part was less variable; we therefore refer to these recombinant genes as Mamu-KIR3DL10A/3DL02 and Mamu-KIR3DL10B/3DL02. The other two recombinant genes consisted of the 5′ end up to exon 3 of Mamu-KIR3DL02 and exon 4 up to the 3′ end of Mamu-KIR3DL08. These recombinant sequences can be distinguished by a synonymous and a nonsynonymous single nucleotide polymorphism (SNP) in the Mamu-KIR3DL08 segment, and will be referred to as Mamu-KIR3DL02/3DL08A and Mamu-KIR3DL02/3DL08B.

The transcripts of the above-mentioned recombinant genes encoded three extracellular domains and a long cytoplasmic tail, suggesting an inhibitory function (Fig. 2C). The recombinant genes were confirmed by independent PCRs and/or segregation into families. All four recombinant genes contained a segment of Mamu-KIR3DL02, suggesting that this gene is highly susceptible to engaging in fusion events. However, recombinant genes containing segments from other inhibitory, activating, or pseudogenes might be discovered when larger populations are studied.

An analysis of a rhesus macaque family comprising 30 animals (Fig. 1) revealed the segregation of KIR gene combinations. Twelve previously unreported Mamu-KIR haplotypes were deduced (Fig. 3A). All haplotypes were confirmed in multiple animals, except for haplotype H12, for which material of informative offspring is missing (sire 95055). This deduced haplotype was, however, confirmed by multiple independent PCRs and, in addition, Sanger sequencing to certify low PacBio read counts for the Mamu-KIR3DL20 and KIR3DS05 genes (<3 PacBio reads; Supplemental Fig. 1A). KIR3DS05 on haplotype H12 could not be defined at the allelic level, but the presence was confirmed on gDNA. Furthermore, haplotype H12 was confirmed in part by a previous study in which the KIR repertoire of sire 95055 was analyzed using Sanger sequencing (22).

The number of KIR genes per haplotype showed remarkable variability, and ranged from 4 to 14 genes, which produced bona fide transcripts. Mamu-KIR3DL20 was transcribed on all 12 haplotypes, whereas Mamu-KIR3DL01 was present on all haplotypes except for haplotype H9. These two genes might be considered framework genes in rhesus macaques, although some previously reported macaque KIR haplotypes defined by conventional characterization methodology seem to lack those genes (22). Eight of the 12 Mamu-KIR haplotypes transcribed a Mamu-KIR2DL04 gene, whereas in humans the orthologous gene is referred to as a framework gene (12). Other common lineage II KIR genes were Mamu-KIR3DL05, Mamu-KIR3DL07, Mamu-KIR3DL08, Mamu-KIR3DL10, Mamu-KIR3DS02, and Mamu-KIR3DS05, which were present in 55–70% of the studied animals. Most of the remaining KIR genes were present in ∼10–30% of the animals, or were absent in this family. The more frequently present genes, like Mamu-KIR3DL01, Mamu-KIR3DL05, Mamu-KIR3DL07, Mamu-KIR3DL20, and KIR3DS02, showed the most allelic variation (Table II), which might indicate a selective pressure on relatively rapidly evolving genes, or that these KIR genes are old entities that accumulated variation over time.

Haplotypes H7, H9, and H11 showed an expansion of the KIR cluster that was characterized by two or three “allelic” copies of a certain gene (Fig. 3A). On haplotypes H7 and H11, Mamu-KIR3DL01 was duplicated, resulting in two allelic Mamu-KIR3DL01 copies that differed at 22 and 12 bp positions, respectively. In view of the number of SNPs between the Mamu-KIR3DL01 copies and the fact that they are located on the same chromosome, designating them as different genes could be considered in the future. Macaque haplotype H9 showed a more extensive expansion, with two allelic copies of Mamu-KIR3DL05 and Mamu-KIR3DS02, and even three copies of Mamu-KIR3DL07. Each allelic Mamu-KIR3DL07 copy varied from the other by at least six SNPs, suggesting that these copies are old entities that accumulated variation over time, and are not the result of recent duplications. The duplicated Mamu-KIR3DL05, Mamu-KIR3DL07, and Mamu-KIR3DS02 copies on haplotype H9 were also confirmed by Sanger sequencing to exclude potential in vitro artifacts generated by the PacBio platform (Supplemental Fig. 1B–D, Supplemental Table I).

In contrast to extended haplotypes, macaque haplotypes H4, H8, H10, and H12 seemed to be contracted haplotypes, with only four, seven, or eight KIR genes present, including a recombinant gene generated by an in-frame fusion event. All of these contracted haplotypes also contained a Mamu-KIR3DL01 and Mamu-KIR3DL20 allele along with a few additional lineage I and/or II KIR alleles. The presence of a fusion gene on all contracted haplotypes might be indicative of recombination events, such as unequal crossing-over and intrachromosomal recombination (Fig. 2A, 2B), which might have caused the deletion of genes that were present on the original KIR haplotypes. At least one copy of a recombinant gene was found in 40% of the animals (Supplemental Table II), which might indicate positive selection for contracted haplotypes containing recombinant genes.

The recombination events in macaques may be indicative of the physical position of KIR genes and, in combination with the previously sequenced macaque KIR haplotype by Sambrook et al. (27), the physical locations of the KIR genes were predicted and illustrated in Fig. 3A. For example, haplotypes H8 and H12 contain a recombinant gene with a 5′ segment of Mamu-KIR3DL10 and a 3′ segment of Mamu-KIR3DL02, suggesting that the latter is located downstream of Mamu-KIR3DL10. The recombinant genes on haplotypes H4 and H10 suggest that Mamu-KIR3DL02 should be localized in front of Mamu-KIR3DL08. Furthermore, Mamu-KIR3DL05, Mamu-KIR3DL07, and Mamu-KIR3DS02 seemed to have been introduced on macaque haplotype H9 as a single entity, indicating that these genes were located next to each other, or are at least in close proximity. Nevertheless, to elucidate the precise KIR gene positions, additional genomic haplotype sequencing and phasing needs to be performed.

To confirm that no KIR genes were missed by PCR amplification at the transcription level, the absence or presence of several frequent KIR genes was also determined at the gDNA level (Supplemental Table II). This approach suggested that no frequent KIR genes had been missed by the three primer sets used for amplification on cDNA. In addition, although typing for the presence or absence of these genes at the gDNA level might be less informative, in a few cases it provided extra information that was needed to assign them to a certain haplotype.

In humans, one generic primer set, mapping to the UTRs, was able to amplify most KIR2D/3D genes, except for KIR3DL3, KIR2DL4, and KIR2DL5. For KIR2DL4, we designed an additional specific primer set to facilitate analysis of this gene at the transcriptional level (Table I). UTR-specific primers for KIR3DL3 and KIR2DL5 could not be designed, and therefore these genes were only analyzed for the presence or absence at the gDNA level. It is known from the literature that KIR3DL3 is a framework gene, and should be present on all KIR haplotypes, whereas for KIR2DL5, two, one, or no copies can be present only on group B haplotypes (12, 45).

A human family consisting of 15 members was selected, which allowed segregation analysis (Fig. 4). The family comprised three generations, and had been founded by five individuals. Approximately 7900 KIR2D/3D reads and 7350 KIR2DL4 reads per individual were obtained, and 9–22,5% of the reads mapped 100% to the human KIR2D/3D allele library containing 907 reported KIR alleles (14). The remaining reads contained single nucleotide gaps, or were partial sequences introduced by PacBio sequencing. No unreported human KIR alleles were discovered, although one known partial sequence was extended (KIR3DL2*011:01, LT934502). Furthermore, a recently reported fusion gene was confirmed, consisting of the 5′ end up to exon 6 of KIR2DL1 and exon 7 to the 3′ end of KIR3DL2 (Fig. 2D, LT963640) (42). The fusion transcript encoded a D1-D2 extracellular segment and a long intercellular tail, suggesting an inhibitory KIR2D receptor. Standard KIR typing kits, which are commonly used for KIR characterization in humans, readily miss recombinant KIR genes as they only type for the presence or absence of gene segments at the gDNA level.

Ten KIR haplotypes from 15 related human individuals were thoroughly defined based on full-length cDNA transcripts (Fig. 3B). Each haplotype encoded from 3 to 10 KIR gene transcripts, including representatives of the framework genes KIR3DL3, KIR2DL4, and KIR3DL2, except for haplotype H4, which lacked a copy of the KIR2DL4 gene. The fourth human framework gene, KIR3DP1, is a pseudogene, and is therefore not amplified at the transcription level.

In this paper, we followed the human KIR haplotype nomenclature conventions of Pyo et al. (46), which was later adapted by Vierra-Green et al. (47). Human haplotypes H1, H2, H5, and H7–H10 represent so-called nonvariable A haplotypes (cA01|tA01), characterized by six KIR genes (on cDNA level), including the framework genes KIR2DL4 and KIR3DL2, and the activating receptor KIR2DS4. Haplotype H8 was identical to the previously reported cA01:009|tA01:010 haplotype, whereas the telomeric region of H7 was identical to tA01:011 (46). Group B haplotypes were represented by haplotypes H3 and H6 (Fig. 3B), containing up to 10 KIR genes, including multiple activating receptors. The telomeric region of haplotype H6 was identical to the previously reported tB01:001 region, and was combined with a cA01 region (cA01|tB01:001). KIR haplotype H3 confirmed the previously reported cB01|cA01 haplotype configuration (46). Haplotype H4 contained only three KIR genes—KIR3DL3, KIR2DL3, and a KIR2DL1/KIR3DL2 fusion gene—the latter suggesting a contracted haplotype. The absence of the framework gene KIR2DL4 on haplotype H4 might suggest a deletion in the telomeric region of this haplotype. Recently, Roe et al. (42) published a contracted haplotype on gDNA that appears to be identical to our haplotype H4, including the above-described recombinant gene (cA04).

In most human disease association studies, KIR genes are typed by determination of their presence or absence on gDNA. However, multiple studies demonstrated that health and disease could be linked to certain KIR alleles instead of to KIR genes (48–50). Therefore, high-resolution characterization of the KIR genes might be clinically beneficial, and could improve future KIR disease association studies. To confirm that our transcriptome characterization approach did amplify all KIR genes, their presence or absence was also assayed at the gDNA level using the Olerup SSP KIR typing kit. All 17 KIR genes that were identified as present on gDNA were also detected at the transcription level, except for KIR2DL5, KIR3DL3, and the pseudogenes (Supplemental Table III). The pseudogenes KIR3DP1 and KIR2DP1 were present in all genotyped individuals. KIR3DP1 is a framework gene and was suggested to be present on all haplotypes, except for haplotype H4, which is assumed to lack the KIR3DP1 gene because of a deletion in the telomeric region as described by Roe et al. (42) (Fig. 3B). KIR2DP1 has been described to be present on cA01 and cB01 regions (46), suggesting the presence of this pseudogene on all haplotypes of the studied human family. The framework gene KIR3DL3 and the group B haplotype-specific KIR2DL5 were found present at the gDNA level in all individuals and on both group B haplotypes H3 and H6, respectively (Fig. 3B; gray boxes). These findings supported the assumption that our transcriptome protocol amplified all KIR genes, except for KIR2DL5, KIR3DL3, and the pseudogenes, and might be beneficial for future KIR disease association studies. In addition, the identification of multiple human KIR haplotype regions that had been reported by others further validated our protocol.

More recently, sequencing and characterization studies provided insights into the complexity of the KIR gene system. However, because of large gene CNV and the high similarity of the KIR genes, conventional sequencing methods hampered the accurate characterization of the complete KIR system at the transcription level. In this study, we describe a comprehensive and relatively fast SMRT sequencing protocol using a PacBio Sequel platform to completely characterize the KIR transcriptomes in human and rhesus macaque families. The power of this approach is demonstrated by the fact that, in a relatively short time, novel KIR genes, alleles, and complex KIR haplotypes were defined by segregation studies in a family setup. A relatively high number of human and macaque recombinant KIR genes were discovered, and seemed to be the result of several independent fusion events. This study allowed comparison of the human and rhesus macaque KIR transcriptomes. Eventually, this may result in a better understanding and interpretation of KIR disease association studies.

A comparison of 10 human and 12 rhesus macaque KIR haplotypes illustrates that both species share highly similar gene systems. Although there are subtle differences, such as different receptor lineage expansion and haplotype organization, both species show extensive allelic polymorphism and gene CNV in their KIR repertoire. Hence, rhesus macaques may provide relevant models to study the impact of KIR genes on health and disease. Human KIR allelic polymorphism seems to have already been broadly mapped, as all the alleles we recovered are documented in a database containing over 900 alleles extracted from numerous population studies (14). In contrast, only 218 Mamu-KIR alleles were reported, but considering that almost half of the total Mamu-KIR alleles discovered in this study were unreported, it is reasonable to suggest that this number is only the tip of the iceberg. Therefore, the extent of allelic KIR polymorphism seems to be at least comparable in humans and rhesus macaques.

KIR gene CNV is observed on the human and rhesus macaque haplotypes, which contained from 4 to 12 or 4 to 14 KIR (pseudo)genes, respectively. In humans, the nonvariable group A haplotypes rarely show CNV, whereas the group B haplotypes can have a variable number of—mainly activating—KIR genes, in part caused by duplications or chromosomal rearrangements (Fig. 3B). In macaques, each haplotype can contain a different number of inhibitory and activating KIR genes, which can be magnified by deletions or insertions as a result of duplications and chromosomal rearrangements (Fig. 3A). In comparison with the human situation, the overall CNV seems to be more extensive in the macaque KIR gene system, and might be explained by coevolution with the expanded MHC class I repertoire in macaques, and by the absence of a haplotypic organization as is observed in humans.

One human and four macaque haplotypes showed signs of contraction. In both species, the contraction of haplotypes was marked by the presence of recombinant genes and the apparent deletion of a haplotype segment. The generation of these contracted haplotypes is most likely mediated by repetitive elements present in the KIR introns (43). For some human recombinant haplotypes, these sequence repeats are identified and characterized as breakpoints that may facilitate chromosomal rearrangements (Fig. 2A, 2B). In macaques, however, these repetitive elements are not yet characterized, but considering the observation of short Mamu-KIR haplotypes, in addition to the presence of recombinant genes, it is likely that the same mechanisms are responsible for KIR gene expansion and haplotype contraction as is observed in humans. The rapid loss and gain of KIR genes, driven by the repetitive sequence elements in the introns, might be an advantageous evolutionary strategy to expand gene variability using an existing gene repertoire, and thereby enhancing pathogen evasion. Hence, recombinant genes are composed of different heads (ligand interaction) and tails (signaling function), which may facilitate the exchange of functionalities and ligand interactions between receptors (Fig. 2C, 2D).

In macaques, in addition to contraction, expansion of KIR haplotypes was also observed. Macaque haplotypes H7 and H11 seem to have expanded by gene duplication, whereas haplotype H9 showed evidence of unequal crossing-over events. However, it is arguable whether the duplication on macaque haplotypes H7 and H11 should be considered as copies of a Mamu-KIR3DL01 gene, because the allelic copies vary at 22 and 12 bp positions, respectively. A sensible nomenclature system for Mamu-KIR3DL01 sequences should be considered. A more extreme expansion is observed on macaque haplotype H9, which contains two allelic copies of Mamu-KIR3DL05 and Mamu-KIR3DS02 and three Mamu-KIR3DL07 copies. Although it is possible that this expansion is explained by multiple gene duplications, it seems more likely that these three duplicated genes are introduced as one entity or tandem, on which Mamu-KIR3DL07 might already have been duplicated.

Although no expanded human KIR haplotypes were found in the present human family studied, other researchers reported KIR haplotypes with gene insertions, similar to the extended macaque haplotype H9 (42, 51, 52).

The expansion and contraction of KIR haplotypes might be a balancing selection for fighting infections on the one hand and for reproductive success on the other. A similar reproductive/immunological trade-off is illustrated by the haplotypic organization in human KIR (12). In humans, expansion and contraction of the KIR region are only observed on group B haplotypes, which indicate that structurally diversifying the nonvariable group A haplotypes is not beneficial. The group B haplotypes, which can contain a variable number of KIR genes, show chromosomal rearrangements, whereby KIR genes are introduced and/or deleted. This diversifying selection might be associated with increasing immune response variability, successful reproduction, loss of unfavorable genes, or the generation of novel fusion genes. However, chromosomal recombination might also be driven by the specific content of group B haplotypes. In this case, sequence elements that are present on group B haplotype-specific genes, such as activating KIR and KIR2DL5, might drive recombination events, without necessarily requiring selective pressure. In macaques, as well as in all other nonhuman primates, there is no haplotypic organization that divides variable and nonvariable haplotype content. All macaque haplotypes seem to be prone to diversifying selection, as great variability in gene content is observed, but association with reproductive success or pathogen evasion has not yet been demonstrated. On an individual level, contracted haplotypes should provide the essential functions of the KIR gene system. Mamu-KIR3DL01 and Mamu-KIR3DL20 are two highly polymorphic and frequently expressed genes, and all four contracted Mamu-KIR haplotypes identified in this study have expression of these genes in common. Three of the studied macaques (R12035, R13084, R14072) are homozygous for contracted haplotypes (Figs. 1, 3A), and do not show signs of an impaired immune system. This might indicate that Mamu-KIR3DL01 and Mamu-KIR3DL20 expression, in combination with a recombinant gene, is sufficient to provide functional NK cell activity.

Previously, two population studies (51, 52) and two smaller studies (42, 53) reported contracted and expanded KIR haplotypes at the gDNA level in humans. In three of these studies, almost twice as many contracted haplotypes were observed in comparison with expanded haplotypes. However, the prevalence of contracted and expanded haplotypes in these populations with European ancestry was only 5–10%. In the selected rhesus macaque family, 4 of the 12 Mamu-KIR haplotypes were contracted, whereas only one showed extensive expansion. In 46% of the animals of the studied rhesus macaque family, a contracted or expanded haplotype was observed. Although a haplotype analysis of a larger rhesus macaque population is required to compare the prevalence of contraction and expansion to that observed in humans, this study suggests that chromosomal rearrangements in the KIR cluster is more common in rhesus macaques than it is in humans. Furthermore, the results indicate that contraction, accompanied by the generation of novel recombinant genes, seems to be more beneficial, or at least has a higher occurrence, than expansion of KIR haplotypes.

In the past, different typing strategies were reported to characterize the KIR cluster in humans and rhesus macaques, including cloning and Sanger sequencing, Roche/454 pyrosequencing, and microsatellite analysis (22–24, 27, 33, 38). However, these methods were often insufficient to assemble full-length allele sequences, as well as being time-consuming and mainly focused on the presence or absence of genes. More recent reports described high-resolution KIR characterization on gDNA by next-generation sequencing on a PacBio RS II platform, a MiSeq platform, or by exome capturing (42, 54, 55). These studies were able to define allelic-level genotypes and to identify novel alleles, and less frequently defined complete or partial haplotypes, and discovered recombinant genes at the gDNA level. Another recent study used SMRT sequencing on a PacBio RS II platform to characterize KIR transcription in an inbred Mauritian cynomolgus macaque population, which is restricted in its genetic diversity, and was able to identify nine novel alleles and to define KIR haplotypes (56). Our studies show that a similar approach is also applicable to outbred human and macaque populations.

Most current clinical methods to characterize the KIR gene system are primarily based on determining the presence or absence of the known human KIR gene segments, and they might miss substantial information such as allele-level typing, CNV, expanded haplotypes, and recombinant genes. Multiple disease association studies, however, illustrated the importance of distinguishing between alleles of KIR genes (48–50). Therefore, a high-resolution KIR characterization approach might be beneficial for future health and disease studies.

In this paper, a method is described to thoroughly characterize the KIR transcriptomes in humans and rhesus macaques, using a relatively fast high-resolution SMRT sequencing protocol on a PacBio Sequel platform. Novel alleles and recombinant genes were discovered, and transcribed haplotypes were defined based on transcription profiles in concert with segregation studies in families. Although sequencing at the transcription level might have minor drawbacks, such as the lack of intron and pseudogene information, it does have serious advantages over sequencing on gDNA; for example, transcriptional modifications can be observed, including splicing of transcripts or intron insertions. These transcriptional modifications might have an effect on the function of the receptors; for instance, when KIR3D transcripts can be spliced to generate KIR2D transcripts, or when transmembrane regions are spliced out, which may result in soluble receptors. In addition, although the PacBio platform does not yet provide quantitative analysis, the number of identified transcripts might be indicative of the expression level of a certain allele. However, this indication might be affected by other factors, such as variegated KIR expression, the MHC class I gene repertoire, and previous pathogen exposure.

In conclusion, this comprehensive sequencing approach can eventually contribute to better understanding and characterizing the KIR gene cluster in different species, thereby improving not only the interpretation of disease association studies but also transplantation and reproduction biology.

We thank D. Devine for editing the manuscript and F. van Hassel for preparing the figures.

The authors have no financial conflicts of interest.