The mammalian immune system provides critical protection against a broad variety of insults but when dysregulated can itself give rise to pathologies. Its proper function and regulation depend critically on a balance of activating and inhibitory signals, which are received and coordinated via cell-surface immunoreceptors and their associated downstream signaling pathways. In vertebrates, the genes underlying these receptors are often found in clusters containing evolutionarily and/or functionally related genes and gene families (1), and these regions are frequently preserved across phylogeny (2).

Although significant attention continues to be paid to the genetic mechanisms underlying the adaptive immune system, less attention has been given to the genetics of the innate immune system, which is essential to the survival of an organism and ultimately a species. Marked differences in the cell repertoire and the functioning of the innate immune system have been documented between closely related species (3, 4), and it has been shown that species within a family can vary in the degree to which they rely on the adaptive vs. innate immune systems to achieve immunocompetence (5). Among the cells with functions related to innate immunity, Natural Killer (NK) cells are notable for their capacity to lyse target host cells, an ability shared only by CD8+ T cells. This capability is significant not only for the elimination of intracellular pathogens but also for the destruction of cancer cells and can be harnessed for therapeutic purposes (6, 7).

Two genomic regions contain clusters of NK cell receptors (NKR) in mammals: the Natural Killer Complex (NKC) and the Leukocyte Receptor Complex (LRC). The NKR genes found in these complexes are structurally related within each complex but differ between the complexes. The NKC contains the killer-cell lectin-like receptors (KLRs) which possess C-type lectin-like extracellular domains, while the LRC genes encode receptors with extracellular ligand-binding domains belonging to the immunoglobulin (Ig) superfamily, such as the killer-cell immunoglobulin-like receptors (KIRs) and the leukocyte immunoglobulin-like receptors (LILRs) (1). Despite their structural differences, representatives of these gene families can bind polymorphic major histocompatibility complex (MHC) class I molecules as ligand. Moreover, they share similar transmembrane and intracellular domains, and therefore can fulfil the same functions in terms of downstream signaling and mediation of NK cell responses (8, 9). Specifically, these gene families code for both activating and inhibitory receptors, and ultimately the balance of activating and inhibitory signals received by a given NK cell determines its activation status and allows the cell to distinguish between self and non-self or altered self (10).

In contrast to the MHC, which has a relatively conserved genomic organization in mammals (11), the genomic regions containing the NKR gene families are evolutionarily flexible and rapidly evolving. In an interesting example of convergent evolution in mammals, different NKR gene families have expanded and diversified in different mammalian species (9, 12). NK cell MHC-I recognition is mediated by receptors encoded in the LRC in humans and higher primates (13), contrary to receptors encoded in the NKC in rodents (14), equids (15), and prosimians (16). The marked difference between mice and rats demonstrates that even within closely related species the gene content of the NKC and LRC can vary significantly.

Within the LRC, the KIR family has been characterized across a variety of placental mammals, with expansion and diversification of the KIR described in some species of primates (platyrrhines and catarrhines) and some artiodactyls (cattle) (9). The carnivores studied to date include domestic dogs, domestic cats, and several pinniped species (17). No expansion of the KIR gene family has been identified in any of these species: the domestic dog genome lacks KIR entirely, while a single KIR gene or pseudogene has been identified in domestic cats, three species of seals, and sea lions.

Less is known about the LILRs, which have mainly been studied in primates (18) and mice (as PIRs, the murine orthologues of LILRs) (19), as well as in several species of Cetartiodactyla, including pigs (20), goats (21), and cattle (22). Like other immunoreceptor gene families, the LILRs comprise both activating and inhibitory receptors (23). Activating LILRAs possess a short cytoplasmic tail and a positively charged amino acid residue in the transmembrane region which associates with immunoreceptor tyrosine-based activation motif (ITAM) containing proteins for signal transduction. Inhibitory LILRBs have long cytoplasmic tails containing 2 to 4 immunoreceptor tyrosine-based inhibitory motifs (ITIMs) (24). In addition, several LILRAB genes with characteristics typical of both activating and inhibitory LILRs have been described in primates (18), and one human LILR without transmembrane and signaling domains is known to be expressed as a secreted protein (25). Both LILRA and LILRB gene products comprise 2 to 6 extracellular Ig domains, depending on the species: either 2 or 4 Ig domains are typical for human LILRs (24), while 6 Ig domains are typical for murine PIRs (26) and can also be found in some bovine LILRs (22). Additionally, group 1 and group 2 LILRs have been defined for pigs, cattle, and goats based on phylogenetic analysis, with human LILRs clustering with the Cetartiodactyla group 2 LILRs (20).

LILRs are expressed in a variety of immune cells, including monocytes, macrophages, B and T lymphocytes, granulocytes, NK cells, mast cells, and dendritic cells (18). Although they have been shown to bind a variety of ligands, their interaction with MHC-I is the best described (27). Other ligands, especially for activating LILRs, include pathogen-associated proteins and host immunomodulatory molecules (28). However, at this time the function of various LILRs is poorly characterized, especially in non-human species. In humans and to a lesser extent in mice, LILRs (and their PIR orthologues) have been shown to play roles in autoimmune and inflammatory diseases, viral responses, neurodegenerative disorders, infectious diseases, and cancer (18, 29). Therefore, they are considered potentially useful as diagnostic markers and as a target for immunotherapies (30–33).

The family Felidae has the potential to be an informative model for comparative immunogenetic studies. Its members include species living in a variety of habitats and expressing a huge diversity of feeding and social behaviors, resulting in exposure to different pathogen pressures. Despite this potential, our knowledge of the immunogenomes of felids is limited, although more attention has been paid to domestic cats than wild felids. Nonetheless, differences in immune function between Felidae species have been documented, such as the increased reliance on innate immune mechanisms in cheetahs compared to a sympatric leopard population (5). Therefore, inferences from one species to another may be inaccurate, and the study of the immunogenome of individual species is warranted. A review of the MHC in Felidae was recently published, contributing to our knowledge of the genes underlying the adaptive immune system in these species (34). However, the genes underlying innate immune receptors, and particularly the LILRs, remain largely undescribed. The current study, therefore, aims to characterize the LRC and particularly the LILR gene family in the felid long-read, chromosome-level genome assemblies currently available.

Until recently, comparative studies of the LILRs have been limited in part by the technical difficulty of sequencing the LRC, which contains a number of closely related genes and is highly repetitive. Such regions are generally difficult to assemble from short-read sequencing data (38, 39), and therefore the LRC is often fragmented or misassembled in short-read assemblies. In scaffold-level assemblies, breaks are especially common within the LILR clusters, and portions of these genes may also be found on small unplaced scaffolds. Long-read assemblies are better able to map such technically challenging regions without breaks; the impact of this method has been such that it was named “Method of the Year 2022” by Nature Methods (40). For these reasons, this study focused on long-read chromosome-level assemblies. However, basecalling errors are more frequent in long-read assemblies (41), despite reported improvements in basecalling for both SMLR (42) and ONT sequencing (43). Thus, some genes with only minor changes in their sequence (1-2 bp indels/mutations) may be erroneously categorized as pseudogenes and in fact may be functional.

The LRC is highly conserved across the ten felid assemblies meeting the inclusion criteria ( Supplementary Figure 1 ) and it is identical in another F. catus genome assembly (GCA_016509815.2). Moreover, no inconsistencies in the content of the LILR sub-region ( Supplementary Figure 2 ) were identified in the short-read or scaffold-level Felidae assemblies ( Supplementary Figure 3 ), nor in the short-read assemblies of other Feliformia species (data not shown).

Within the Carnivora more largely, chromosome-level long-read assemblies were available for several species of both the Canidae and Mustelidae families, allowing the intrafamily variability of the LRC to be assessed and compared with the Felidae. Because these assemblies met the same inclusion criteria as the focal felid assemblies, concerns that apparent variability may be the result of different sequencing technologies are reduced. In addition to these two families, a long-read chromosome-level assembly was available for the California sea lion; this is the only other assembly of a Carnivora species meeting these criteria, and it has been included to broaden the view of the order. The overall architecture of the LRC was again consistent across species, suggesting that the conserved structure of the LRC observed in the Felidae is typical for the order Carnivora. However, two partial exceptions should be noted: the wolf and the American mink.

In the wolf, the block of genes spanning from FCAR to LILRB3 appears to be duplicated. Although the sequences of this region and its duplicate are very similar, they are not identical. The cognate pairs of the (pseudo)genes often contain 2 to 3 substitutions per gene and differ in length by up to 52 base pairs. LILRB3 and LILRB3⁺ , the only pair of putatively functional LILRs resulting from the duplication, differ sufficiently in their sequences that on phylogenetic analysis LILRB3⁺ clustered more closely with the orthologous gene in the dingo and the dog than with wolf LILRB3. However, whether this duplication is a true feature of the wolf genome or an artifact of the assembly remains to be determined.

The apparent division of the LRC in the American mink into two subregions separated by approximately 15 megabases represents a significant divergence. Previous research on the phylogeny of the LRC has not, to the best of the authors’ knowledge, described such a division in any other mammalian species (2). Nonetheless, the fusion of two chromosomes which are separate in many Mustelidae karyotypes (including European mink, forest polecat, lesser weasel, mountain weasel, Japanese sable, striped polecat, and European badger) into the single American mink chromosome 7 has been elucidated by G-banding and FISH data (44). The genes present on the respective chromosomes in the examined assemblies confirm the fusion of M. meles chromosomes 8 and 19 into N. vison chromosome 7. Graphodatsky et al. (44) propose an ancestral Carnivora karyotype (ACK), with ACK chromosome 19 corresponding to chromosome E2 in the Felidae and to the p arm of N. vison chromosome 7, which corresponds with the locations of the LRC in these species, and confirming the homology described in an earlier FISH study (45). Thus, it may be that during the karyotype evolution that produced N. vison chromosome 7, chromosome breakage and reassembly split the LRC into two regions. Unfortunately, no other assemblies of Neogale spp. are available for comparison; the only other assembly available for the genus is a second N. vison genome, which is assembled only to the scaffold level and is insufficient to confirm or disprove the division of the LRC. The closest relatives for which assemblies are available are the Mustela spp.; however, the chromosomal rearrangement affecting N. vison chromosome 7 did not occur in the Mustelae (44). This two-part arrangement of the American mink LRC could allow for the generation of additional diversity in the LILRs as they would be less likely to be inherited en bloc. In this context, it is interesting that the mink possesses the largest repertoire of LILR (pseudo)genes of any of the studied Mustelidae. This possible difference in the genomic organization of the LRC between closely related species demonstrates the importance of studies conducted at the species level in rapidly evolving regions such as the LRC.

The SMLR felid assemblies each contain seven orthologous LILRs. Although we hypothesized that there could be differences in the LILR content between the Felinae and Pantherinae, no such differences were identified. The preserved localization of gene components (introns, exons, UTRs) in the VISTA graphic further attests to the conserved nature of the LILR family in Felidae ( Supplementary Figure 2 ). In the case of the Canada lynx, pseudogenes should be considered in light of the identified anomalies in the assembly. Similarly, the clouded leopard assembly was recently made available and is not yet expected to be in its final form. In the tiger, some pseudogenes again may be an artifact or the result of interindividual variability, as LILRA3 is putatively functional in a recently released assembly based on linked short reads and using Hi-C scaffolding (GCA_024034525.1) (46). The individual used for the SMLR reference assembly was a cross-species hybrid (tiger X lion), and it is possible that the challenge of separating the two sequences led to some misassembly in technically challenging regions like the LRC. The cheetah, on the other hand, has the most pseudogenes (three), and interestingly the cheetah immune system is known to have several specific features. The variation in the cheetah MHC has traditionally been considered to be notably low (47). While much higher genetic diversity within MHC I (48) and MHC II (49) loci has been described recently, cheetahs still show MHC II-DRB diversities lower than other large felids like Bengal tigers (50), Eurasian lynx (51), and leopards (52). Cheetahs seem to achieve immunocompetence through higher constitutive innate immunity but lower induced innate and adaptive immunity compared to a sympatric leopard population (5). It is tempting to hypothesize that the higher number of LILR pseudogenes in this species could be related to these specificities, but focused research of LILR expression and function is needed to draw any conclusions.

Differences in the LILR gene content of the Felid assemblies not meeting the inclusion criteria may to some extent reflect the sequencing technology. In the bobcat and the leopard, the two SMLR scaffold-level assemblies, orthologues of the seven described LILRs were identified. This was also the case for the snow leopard assembly, which was generated by linked short reads. The manul assembly is also a scaffold-level long-read assembly but based on ONT rather than SMLR, and it also contains these seven orthologues, although only two appear to be functional. The cougar assembly, generated by a combination of short-reads and ONT, contains many LILR fragments but only two full genes, while the two short-read assemblies with the LRC contained on a single scaffold, the black-footed cat and the jaguar, each contain two LILRs and no fragments within the LRC. However, the number of Ig domains identified in a BLAST of the WGS in both cases far exceeded the number of such domains seen in the LRC along with the X chromosome, and in both cases was similar to the number of Ig domains found by BLAST of the WGS in the domestic cat. This suggests that in these assemblies some LILRs may be found on small unplaced scaffolds. This correlation of the number of identified LILRs with the use of different sequencing technologies suggests that long-reads are indeed more successful at resolving the LILR subregions.

Regarding the Carnivora more broadly, the overall LILR content is similar between the Felidae (seven LILRs), the Canidae (six to eight), and Z. californianus (nine). A possible exception is the wolf, which may have a total of 11 LILRs if the duplication described is a true feature vs. eight if it is taken as an artifact of the assembly. The Mustelidae possess more LILR genes, ranging from 12 (Eurasian otter) to 23 (American mink) in the species studied. The number of LILR fragments is also higher in the Mustelidae. However, this does not translate into a similarly higher number of functional LILRs in the mustelids: on average they possess 6.75 putatively functional LILRs, compared to 6.3 in the felids. Indeed, Felidae possess a notably higher ratio of functional LILRs to pseudogenes (90.00%) compared to the Mustelidae (41.54%) or the Canidae (57.50%; 59.46% with exclusion of the wolf duplication). Due to the so far poor characterization of LILR functions, it is difficult to make a link between interspecific variation in the LRC gene count and potential selective pressures leading to such a differentiation.

The intrafamily variability in the LILR complement was lowest in the Felidae. The other end of the spectrum is again represented by the Mustelidae, which show the largest variability between species in terms of the total number of LILR loci, the total number of putatively functional genes, and the ratio of putatively functional to pseudogenized LILRs.

The significance of the enlarged family of LILRs and its interspecies variability in the mustelids remains to be elucidated. One hypothesis may be that in the Mustelidae the LILRs fulfill the function of an expanded NK cell receptor family binding MHC class I ligands, similarly to the KIR in humans (53) and KLRA in mice (54). In this interpretation, many genes and pseudogenes could be the result of changing pathogen pressures leading to the rapid evolution and inactivation of LILR genes. However, the LILRs bind the α3 and β2 microglobulin domains of MHC-I rather than the highly variable α1 and α2 domains where KIR and KLRA bind (24), and this may limit their capacity to fulfil a similar function, although they are known to be capable of distinguishing between self and non-self MHC-I in mice (29, 55). Alternatively, the expansion of the Mustelidae LILR family could be more a feature of the evolution of the chromosomal region than the result of selective pressures on the LILRs, and pseudogenes may have arisen as nonfunctional sequences through mechanisms such as the duplication of neighboring sequences.

A putatively functional novel Ig-like gene was found between GP6 and RDH13 in all studied felids except the Canada lynx (caveat again, some misassembly in this part of the Lynx assembly is suspected), as well as in the short-read or scaffold-level Felidae assemblies. A pseudogene or fragment homologous and syntenic to this gene was identified in all studied Feliformia and Canidae species. Comparably, only a signal peptide and Ig-like group was found at this locus in the Californian sea lion. No Ig-like gene was found here in any of the Mustelidae species. Thus, the Felidae are so far the only Carnivora family identified carrying a putatively functional copy of this gene.

Phylogenetic analysis suggests that this gene is related to the novel Ig-like genes recently identified in pigs and subsequently in other species including cattle, sheep, goats, horses (20) and camels (56). The gene content of this Ig-like family varies from one (camels) to seven (cattle) genes. In the previously described cases, these genes are located between NLRP2 and NLRP7, which differs slightly from the locus identified in Felidae. No genes or immunoglobulin domain coding sequences were found in any felid between NLRP2 and NLRP7 in our study.

Previously published studies characterizing the LRC and the LILRs for multiple species within a family are limited. Those which can serve for comparison to the Felidae are the Hominidae, represented by humans, chimpanzees, bonobos, gorillas, and orangutans (18); the Bovidae, represented by cattle (Bos taurus) (22, 57) and goats (Capra hircus) (21); and the Camelidae, represented by the Bactrian camel (Camelus bactrianus), the dromedary (Camelus dromedarius), and the wild camel (Camelus ferus) (56).

The complement of genes in the LRC in each of these families is similar to that described for the Felidae, but the relative orientation of some blocks of framing genes differs in the Bovidae, and in some places differs even between cattle and goats (21, 22). The architecture of the LRC in the Old World camels is most similar to that of goats, with the block of genes from LAIR1 to CDC42EP5 inverted in comparison with the Carnivora. In all cases the LILRs remain divided into two clusters.

The LILR complement in the Hominidae and the Bovidae appears to be larger and more variable than in the Felidae. The Hominidae possess 13 (human, gorilla, orangutan) or 14 (chimpanzee, bonobo) LILRs each, of which seven (gorilla) to 12 (chimpanzee) are putatively functional (18). Relative to the Felidae, this is a larger number of both total and putatively functional LILRs, as well as greater intrafamily diversity in the number of LILRs and the proportion of functional genes. In the Bovidae, sixteen LILRs and the closely related FCG2R gene have been described in cattle (22), comprising five activating and seven inhibitory receptors (70.6% potentially functional genes) while five soluble forms may represent pseudogenes. Eight LILRs and FCG2R have been described in goats, of which four are activating receptors and two are inhibitory (66.7% putatively functional genes) (21). Comparing just these two species shows higher variability in the LILR complement of the Bovidae than was found in the Felidae.

The Camelidae LILR family is comparable to that of the Felidae, with six (Bactrian camel, dromedary) or seven (wild camel) LILRs present in each species. These are all putatively functional in the SMLR assembly of the wild camel, while only three (C. dromedarius) or four (C. bactrianus) are functional based on targeted resequencing in the domestic species (56). This is one of the most similar families to the Felidae in terms of the LILR gene content, and it is notable that more closely related families both to the Felidae (e.g., the Mustelidae) and to the Camelidae (e.g., the Bovidae) are more divergent in their complement of LILRs. Thus, the analysis of an evolutionarily close species is not a reliable predictor of the LILR complement in an unstudied species.

Importantly, the methodology employed in the identification of LILRs may impact study findings. A comparison of the LILR gene content across 48 mammalian species was carried out using OrthoFinder augmented by a reciprocal best BLASTn hit search (3). A smaller complement of LILRs was reported for each of the three Carnivora species in that study than is presented here. Discrepancies also exist between the Hilton study and other studies discussed above, including differences in the number of LILRs identified in each of the higher primates (18) and in goats (21). In most cases, Hilton et al. counted, based on the automated annotations, fewer LILRs than other studies of the same species. These annotations are often inaccurate for tandemly duplicated gene families. In a parallel to the identification of fewer LILRs in the short-read Felidae assemblies in this study, it may be that short-read assembly quality in the region of the LRC inhibited the identification of LILRs in some cases. In species for which a long-read assembly is not available, a BLAST search of the WGS followed by targeted resequencing could be carried out to augment the findings from short-read assemblies but may still overlook some genes, e.g., those split onto two contigs. This underscores the importance of assembly quality in assessing the LRC and particularly the LILRs, as well as the challenges of comparing between species and/or studies relying on different assembly technologies.

Several further steps may contribute to our understanding of the characteristics of the LILRs in Felidae and other families. Targeted resequencing of LILR genes to confirm their sequences and assess allelic variants would provide information about both the inter- and intraspecies variability of the LILRs and would allow subsequent selection analyses to identify sites under purifying and diversifying selection. In addition, expression studies would confirm the expression status of these genes in different tissue/cell types and may illuminate connections between their specificities in A. jubatus and known features of the cheetah immune system. Further, an expanded study of the LILR gene family in other Carnivora species, once long-read assemblies become available, would give a more complete picture of the order and confirm or repudiate the trends seen in the species sampled here. In light of the known differences in the LILR gene family of closely related species, such studies would be of value.

In conclusion, the genomic architecture of the LRC is highly conserved across the Felidae. This same overall structure is also largely conserved in the other studied carnivore species, with the major exception of N. vison, in which a chromosomal rearrangement appears to have split the LRC into two separate regions. There also appears to be a duplication of a portion of the LRC in C. lupus. The organization of the LRC is similar in the more distant Hominidae, Bovidae, and Camelidae families, with the caveat that one or more blocks of genes appear to be inverted relative to their position in the Carnivora LRC. The LILR gene content is also conserved within the Felidae, with a total of 7 orthologous (pseudo)genes identified in each species, all of which are functional in most studied felids. This is similar to the number of LILRs found in the Canidae (6 to 11) but lower than that found in the Mustelidae (12 to 23), and the percentage of pseudogenes is higher in both of these families. Among the studied carnivores, the variability of the LILR gene content between species of a family is lowest in the Felidae and highest in the Mustelidae. Overall, the phylogenetic tree of LILRs reflected the evolutionary relationships among the species, although there were some minor exceptions within the Felidae. Within the carnivores, presumed LILR orthologues clustered together more closely than with other LILRs from the same species but only within families. Two lineages of LILRs were identified in the Carnivora as previously described in Cetartiodactyla. A gene that clustered with the novel Ig-like gene described by Schwartz and Hammond (20) was also identified but was putatively functional only in the Felidae. In comparison to more distant families, the total number of LILRs is higher in the Hominidae, similar in the Camelidae, and varies significantly between species in the Bovidae. The within-family variability in the LILR complement is larger in the Hominidae and the Bovidae compared to the Felidae and is similar in the Camelidae. This indicates that evolutionary closeness is not a good predictor of similarity in the LILR gene family, potentially due to the action of selective pressures and the evolutionary flexibility of the LRC. As high quality long-read assemblies become available for more species, the characterization of their LRCs will increase our knowledge of the interspecies variability of this region and the LILR gene family. Significant gaps remain in our knowledge of LILRs, and comparative studies may improve our understanding of these important genes in humans as well as informing our understanding of the phylogeny of the immune system and potentially offering targets for clinical diagnostics or treatment in the studied species.

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding author.

AJ: Searched the available genomes, mapped the LRCs, performed bioinformatic analyses, extracted CDS for phylogenetic analysis, and drafted and edited the manuscript. JF: Participated in the study design, cross-checked all sequences, performed phylogenetic analyses, and edited the manuscript. PB: Conceptualization and design of the project and edited the manuscript. PH: Conceptualization and design of the study, organized collaborations, and edited the manuscript. All authors contributed to the article and approved the submitted version.