- 1Department of Surgical and Radiological Sciences, School of Veterinary Medicine, University of California, Davis, Davis, CA, United States
- 2Department of Anatomy, Physiology and Cell Biology, School of Veterinary Medicine, University of California, Davis, Davis, CA, United States
- 3Department of Pathology, Microbiology and Immunology, School of Veterinary Medicine, University of California, Davis, Davis, CA, United States
- 4Department of Mathematics and Statistics, University of Guelph, Guelph, ON, Canada
Introduction: The domestic cat (Felis catus) is a valued companion animal and a model for virally induced cancers and immunodeficiencies. However, species-specific limitations such as a scarcity of immune cell markers constrain our ability to resolve immune cell subsets at sufficient detail. The goal of this study was to characterize circulating feline T cells and other leukocytes based on their transcriptomic landscape and T-cell receptor repertoire using single cell RNA-sequencing.
Methods: Peripheral blood from 4 healthy cats was enriched for T cells by flow cytometry cell sorting using a mouse anti-feline CD5 monoclonal antibody. Libraries for whole transcriptome, αβ T cell receptor transcripts and γδ T cell receptor transcripts were constructed using the 10x Genomics Chromium Next GEM Single Cell 5’ reagent kit and the Chromium Single Cell V(D)J Enrichment Kit with custom reverse primers for the feline orthologs.
Results: Unsupervised clustering of whole transcriptome data revealed 7 major cell populations - T cells, neutrophils, monocytic cells, B cells, plasmacytoid dendritic cells, mast cells and platelets. Sub cluster analysis of T cells resolved naive (CD4+ and CD8+), CD4+ effector T cells, CD8+ cytotoxic T cells and γδ T cells. Cross species analysis revealed a high conservation of T cell subsets along an effector gradient with equitable representation of veterinary species (horse, dog, pig) and humans with the cat. Our V(D)J repertoire analysis identified a subset of CD8+ cytotoxic T cells with skewed TRA and TRB gene usage, conserved TRA and TRB junctional motifs, restricted TRA/TRB pairing and reduced diversity in TRG junctional length. We also identified a public γδ T cell subset with invariant TRD and TRG chains and a CD4+ TEM-like phenotype. Among monocytic cells, we resolved three clusters of classical monocytes with polarization into pro- and anti-inflammatory phenotypes in addition to a cluster of conventional dendritic cells. Lastly, our neutrophil sub clustering revealed a larger mature neutrophil cluster and a smaller exhausted/activated cluster.
Discussion: Our study is the first to characterize subsets of circulating T cells utilizing an integrative approach of single cell RNA-sequencing, V(D)J repertoire analysis and cross species analysis. In addition, we characterize the transcriptome of several myeloid cell subsets and demonstrate immune cell relatedness across different species.
1 Introduction
The domestic cat, Felis catus, is a valued companion animal and the second most popular pet in the United States. As of the American Veterinary Medicine Association (AVMA) pet ownership statistics, there are upwards of 46.5 million households with at least one cat in 2023 with ownership trends rising in the last decade from 25% in 2016 to 29% in 2022 (1, 2). In addition to the growing need for feline veterinary care, the cat is an important model for infectious diseases such as virally-induced cancers (Feline leukemia virus, FeLV) and virally mediated immunodeficiency (Feline immunodeficiency virus, FIV) (3–6). Despite the value of the cat as a companion animal and infectious disease model, there are species specific research limitations including few cat-reactive immunophenotyping reagents and a poor knowledge base of immune cell markers and behavior.
In the advent of deciphering the feline immune system, many studies characterizing feline immune cells have relied strongly on antibody-based assays including flow cytometry and immunohistochemistry. This has allowed the study of immune cells in health and disease conditions such as FeLV or FIV (7–9). This approach depends on the assumption that feline immune cells are similar to other species. However, feline immune mediated diseases behave very differently from those in other small animals such as dogs, exhibiting unique etiologies and pathogenesis (10). In a cross-species context, innate immunity shows relatively high evolutionary conservation; however, this conservation dwindles with progression to adaptive immunity and complex immune phenotypes (11). Additionally, there is growing evidence demonstrating that evolutionary relationships do not translate to immune transcriptional or cell type relationships. Such an example is the mouse being a better model for human immune cells than the macaque (12). Thus, understanding the species-specific system in the context of our growing evolutionary database is desirable and there is a strong need to characterize the heterogeneous peripheral immune cell population in the cat, especially diverse populations such as T cells.
In efforts to investigate feline specific T cell diseases, our laboratory has contributed to the growing feline specific reagent and immune cell knowledge base including T cell receptor (TCR) repertoire analysis (13). On this trajectory, single cell RNA-sequencing is a logical next step. It enables the efficient large-scale capture and analysis of heterogeneous tissues, overcoming the challenges of species-specific reagent assays (14–17). Consequently, there has been a bloom in single cell RNA-sequencing studies of the peripheral leukocytes in a variety of mammalian species including the dog, cow, horse, and pig (18–22). Additionally, single cell data facilitates the study of cross species variations in immune cell type and evolution of these cells (12). ScRNA-seq of feline circulating leukocytes have been previously performed demonstrating the presence of 5 major cell types (T cells, B cells, NK cells, monocytes, dendritic cells) in the context of other non-model species (23). However, in-depth characterization of key subtypes and their marker conservation has yet to be determined. Thus, the goal of this study is to characterize the heterogeneity of circulating feline T cells and other captured leukocytes utilizing CD5 flow cytometry enriched scRNA-seq and V(D)J repertoire analysis in clinically healthy domestic shorthair cats between the ages of 6 months and 9 years. CD5 was chosen as a selective T-cell marker given the lack of availability of anti-feline CD3 antibodies. Additionally, we performed a cross species transcriptomic integration to assess feline T cells in the context of 4 other mammalian species - dog, horse, human, and pig. Our results provide a foundation for feline immune transcriptomics.
2 Results
2.1 Single cell atlas of the feline circulating immune cells
ScRNA-seq data obtained from 4 representative healthy cats of different ages consisted of 30,073 quality cells. These age groups were selected as they capture the major feline life stages: 6-month-old juvenile (6MO), 1-year-old mature (1YO), 4-year-old adult (4YO), and 9-year-old aged (9YO). The average number of cells captured per cat was 7,518 with an average of 3,655 transcripts per cell (Figure 1A; Supplementary Presentation 1). Due to the poor annotation of the cat reference genome, we utilized a custom homologous mapping script based on ENSEMBL homologous mapping between the cat and human to improve gene readability. Our script increased the number of gene symbols mapped from 13,886 to 16,936, improving the quality of cell type assignment and biological interpretation. Following batch correction via reciprocal principal component analysis, unsupervised clustering revealed the presence of 21 clusters (Figure 1B). A majority of cell types are equally represented from each sample (Figures 1C, D). A quality control violin plot for clusters can be found in the Supplementary Materials (Supplementary Presentation 1).
Figure 1. scRNA-seq atlas of CD5+ enriched circulating feline immune cells revealed 7 major types across 21 clusters. (A) Table of cell counts from each age group. (B) UMAP plot demonstrating the unsupervised clustering results for feline circulating immune cells from 4 healthy cats of different age groups. (C) UMAP plot of global clustering colored by age. (D) Table of cell type counts across clusters by age. (E) Dot plot demonstrating cell type specific marker expression of the 7 major cell types. (F–O) Feature UMAP plots demonstrating expression profile of key markers for each of the 7 different cell types.
These clusters were categorized into 7 major immune cell types based on canonical marker genes from multi-species PBMC single cell literature, top marker genes and over-enrichment gene ontology (GO) (Supplementary Tables 1, 2). The most numerous cell type identified was T cells called by high expression of T cell receptor complex genes (CD3D, CD3E) and species conserved T cell markers LCK, ITK and TCF7 (12, 23) (Figures 1E–G; Supplementary Presentation 1). Across most T cell clusters, top identified GO terms of biological processes (BP) included T leukocyte processes of differentiation, migration, activation, and effector functions (Supplementary Table 2).
The second most abundant population was neutrophils marked by neutrophil-derived proteases including elastase (ELANE), granulocyte colony stimulating factor receptor (CSF3R) and calprotectin (S100A8/S100A9) (24–26) (Figures 1E, H). In most species, neutrophils are not captured in the PBMC layer during density centrifugation (27). Interestingly, our single cell data set captured a significant number of neutrophils from all but the 1YO sample (Figure 1D). This is the first description of feline neutrophils at the single cell transcriptomic level and can provide insights into biological processes of these abundant cells.
A substantial number of monocytic cells were captured with marker genes including scavenger receptor for haptoglobin-hemoglobin complexes CD163, antigen presenting cell co-stimulatory molecule CD83, lipopolysaccharide detector protein CD14, phospholipase B domain containing 1 (PLBD1) and colony stimulating factor receptor (CSF1R) (28, 29) (Figures 1E, I, J). A relatively small but heterogeneous population of B cells were captured which were identified by conserved markers including B cell surface molecule MS4A1 and Ig-alpha protein of the B-cell antigen component CD79A (23) (Figures 1E, K, L). This population also expresses immunoglobulin genes (IGHG3, IGHM) indicating a sub-population of plasma cells (12) (Figure 1E).
Unexpectedly, we captured plasmacytoid dendritic cells (pDC) and mast cells, which represented rare but discrete cell populations. pDCs were identified by expression of Fc fragment of gamma immunoglobulin (FCRLA), dendritic marker PLAC8B, and immunoglobulin lambda constant 7 (IGLC7) and lack of expression of B cell canonical markers MS4A1 and CD79A (30–32) (Figures 1E, L, M). pDcs clustered more closely with B cells than other myeloid cells. This is due to close lineage associations and the shared expression of integral markers FCRLA and IGLC7 among others (33) (Figures 1E, L). Currently, there is debate regarding the reclassification of pDCs as innate lymphocytes and this analysis further supports such an argument (34) (Figure 1B). Mast cells were marked by IgE receptor MS4A2 and a key component of mast cell protease, CPA3 (35, 36) (Figures 1E, N).
Lastly, we captured a population of non-leukocytic immune cells which were annotated as platelets due to high expression of platelet specific proteins PPBP, VWF and GP9 (37) (Figures 1E, O). This cluster also exhibited an extremelylow average count of unique RNA and total RNA, which is consistent with a highly specialized cell type lacking a nucleus (Supplementary Presentation 1).
2.2 Feline circulating T cells
To further characterize T cell subsets, we independently analyzed T cell clusters previously described at the global level (clusters 0,2,9,1,6,4,5,14,8,12, and 15). Unsupervised clustering of T cells revealed the presence of seven T cell clusters, which were unrecognizable at the global level (Figure 2A). Largely, these clusters represent greater T cell phenotypes along increasing pseudotime values (proxy of differentiation or lineage) (Figures 2B, C). The subtle transcriptional differences, as noted through overlapping differentially expressed genes, make effective annotation difficult (Supplementary Table 3). Nonetheless, we were able to identify several canonical subsets.
Figure 2. CD5+ enriched T cells segregate into naïve T cell subtypes and effectors. Unsupervised clustering of T cells reveals 7 subtypes. (A) UMAP of the scRNA-seq atlas of T cells. (B) UMAP of T cells colored by Pseudotime. (C) Dot plot of marker genes expressed by T cell type. (D) Table of cell type frequencies across ages by cluster. (E) Feature plots of expression and co-expression of TRG and TRDC. (F–Q) UMAP of T cells colored by classical T marker genes defining CD4/CD8 status, naive (SELL, CCR7), effectorness (ANXA2, LGALS3, GZMK, PRF1), terminal differentiation (CCL5) and T exhaustion (TIGIT, PDCD1).
T cell subsets included CD4+, CD8A+/CD8B+ and γδ T cells and were equally represented in all age groups (Figures 2D–H). A great majority were identified as naive cells, distinguished by their high expression of CCR7 and L-selectin (SELL) (38–40) (Figures 2I, J). These naive populations segregated into two larger groups which could be distinguished by differential expression of CD4 and CD8A/CD8B (Figures 2F–H). The CD8A+/CD8B+ naive population represents a divergence from the traditionally observed distribution of naive T cells. Most naive T cells in mammalian species are CD4+ (12). However, a large population of CD8A+/CD8B+ naive cells has been described in non-model species such as the horse (22). This highlights a potentially influential role of CD8A+/CD8B+ T cells in cats.
The clusters 3 and 4, adjacent to the large naive populations, were identified as effector cells (TEM) (Figure 2A). These clusters show a gradual upregulation of T effector molecules including S100A11, galectins (LGALS1, LGALS3), annexins (ANXA1, ANXA2), and granzymes (GZMA, GZMK) (41, 42) (Figures 2C, K–M). Szabo et al. identified CCL5 as a critical marker for highly differentiated effector T cells across tissues as well as its antagonistic expression to the naive marker SELL (41). In our data, Cluster 6 expressed the highest levels of cytotoxic markers including CCL5 and exhibited the highest pseudotime score with no naive marker expression (Figures 2C, M–O). This indicates that this cluster is most consistent with the annotation of a terminally differentiated T cell cluster (CD8+ cytotoxic). Additionally, some effector cells expressed genes associated with T cell exhaustion such as PDCD1 and TIGIT indicating the presence of a few of these cells in circulation (Figures 2C, P, Q).
Cluster 5 represented γδ T cells as inferred by the overlapping expression of the TCR genes TRG and TRDC (13) (Figure 2E). Within this population, there was an expression gradient of naive as well as effector markers, indicating subtypes within the γδ population (Figures 2I, K, L). Although these T cell clusters represent biologically relevant groups statistically, the literature across various species identifies significant challenges associated with clustering of transcriptionally similar but heterogenous T cell populations (43). The γδ cluster exemplifies this challenge in that the cluster’s defining features overlap between two different but not mutually exclusive phenotypes - γδ recombination and effectorness/naiveness. These cells exhibit overlapping expressions of TRG and TRDC, indicating their TCR recombination-based categorization. An overwhelming majority of these γδ T cells are naive given their expression of naive marker SELL (Figures 2I, J). However, some cells appear to express features of effectorness (Figures 2K–N). This suggests some limitations of transcriptionally driven phenotyping of T cells. Nonetheless, this investigation represents a significant contribution towards deciphering the biology of feline T cells.
2.3 Feline circulating T effector cells
As described previously, subtyping T cells with overlapping phenotypes represents a significant challenge, which was further exacerbated by the low number of effector T cells in our dataset. Hence, we clustered effector T cells (clusters 3,4) independently, revealing 13 effector clusters (Figure 3A). Many of these clusters represent minimally divergent biological states as reflected by the marked overlap of top differentially expressed genes (Figure 3B; Supplementary Table 4). However, direct visualization of known phenotypes demonstrate that the effector populations captured exhibit phenotypes of human TH1, TH2, TH17 and Treg cells (Figures 3C–J) (41, 42, 44). Gene sets for each module calculation are provided in the Supplementary Materials (Supplementary Presentation 1). These phenotypes do overlap and are not consistent with identified clusters. Nevertheless, the higher scores among closely positioned cells likely indicates the presence of these subtypes. However, a greater sample size of effector cells would be necessary to better resolve subtypes and their trajectory.
Figure 3. Feline Effector T cells (TEM) segregate into transcriptionally similar clusters and reveal the presence of helper T phenotypes. Unsupervised clustering of TEM reveals 13 subtypes. (A) UMAP of scRNA-seq atlas of TEM. (B) Dot plot demonstrating up to 3 top differentially expressed genes for each cluster determined via Wilcoxon rank sum testing (Adj P <0.05). (C, E, G, I) UMAP of TEM colored by helper T subtype gene modules. (D, F, H, J) Expression UMAP of a representative gene from each helper T gene module presented in parallel.
2.4 Cross species analysis of circulating T cells
To understand how the captured feline T cell populations compared to those of other domestic animal species and humans, we performed an integrative cross-species analysis. We selected 4 species for which peer-reviewed annotated scRNA-seq data with over 15,000 cells was available. ScRNA-seq data from three other veterinary species (dog, horse and pig) and humans were chosen as the reference for homologous genome mapping (18, 21, 22, 45). In total, we integrated 95,366 cells and identified 13 clusters across the 5 species (Figures 4A–F). These clusters were functionally annotated based on markers for T cell identity, naive phenotype, effector phenotype and cytotoxic phenotype in addition to their differential gene expression analysis (Figures 4G, H; Supplementary Table 5).
Figure 4. Cross-species integrative analysis of T cells reveals missing cytotoxic effectors in the cat. Unsupervised clustering revealed 12 clusters across 5 species. (A) UMAP of scRNA-seq atlas of T cells. (B–F) UMAP split by species- dog, horse, cat, human and pig. (G) UMAP of T cells colored by T cell phenotype. (H) Dot plot of marker gene sets for T cell subtypes. (I) Percentage bar chart of clusters stacked by species. (J) UMAP colored by CD5 expression. (K) Bar chart of percentage CD5+ cells per cluster. (L) Bar chart of average CD5 expression across cells in each cluster. (M) Scatter plot of percentage CD5+ cells per cluster versus number of cat cells within the corresponding cluster.
Although feline cells were overrepresented in some naive clusters such as 3 and 8, there was limited representation in clusters 5, 6, and 9 particularly (Figure 4I). To investigate whether this finding was associated with our enrichment protocol, we assessed CD5 transcript expression across all clusters (Figure 4J). Our results suggest that the number of CD5+ cells and average CD5 expression were lowest in clusters 5 and 9 (Figures 4K, L). Additionally, we saw a direct relationship between the number of feline T cells and the CD5+ cell count in each cluster. Particularly, clusters 5, 6, 9, 11, and 12 were most underrepresented (Figure 4M). These clusters were identified as cytotoxic effector cells (Figure 4H). This highlights a limitation of our CD5 enrichment protocol as it may select against specific effector clusters such as CD8+ cytotoxic cells. Much of our current understanding of CD8+ feline T cells is from flow cytometric studies with no further subtyping due to lack of reagents. Thus, further investigation is necessary to validate whether this low frequency is biological or technical.
2.5 V(D)J recombination
To comprehensively characterize T cell receptor (TR) expression, we sequenced TR alpha (TRA), TR beta (TRB), TR gamma (TRG) and TR delta (TRD) transcripts using two dedicated 5’ VDJ libraries generated by amplifying TR transcripts with universal forward primers and C gene-specific reverse primers. The majority of CD4+ naïve, CD4+ TEM and CD8+ naïve T cells expressed TRA/TRB transcripts only, while 60.6% of CD8+ cytotoxic T cells additionally expressed TRG transcripts (Figure 5A). Compared to other αβ T cell subsets, CD8+ cytotoxic T cells exhibited differential TRAV gene usage, and to a lesser degree, TRBV gene usage (Figure 5B). CD8+ cytotoxic T cells with TRG expression utilized distinct TRA and TRB V/J pairings, that were not used by CD8+ cytotoxic T cells without TRG expression or other αβ T cell subsets (e.g. TRAV23/TRAJ25, TRAV25/TRAJ41, TRBV25/TRBJ2-6) (Figure 5C). TRG transcripts in these cells were skewed towards a junctional length of 16 amino acids, had dominant usage of TRGV2-2/TRGJ2-2 genes (Figure 5C) and converged on the conserved motif ‘CAAWDPRGYGWAHKVF’ (Figure 5D). Additionally, these cells exhibited a significantly reduced combinatorial diversity between TRA and TRB chains, with three dominant pairings (Figure 5E).
Figure 5. (A) TCR chain expression of T cell subsets (top). TRA/TRB transcripts dominate in all αβ T cell subsets except for CD8+ cytotoxic T cells, which primarily express TRA/TRB/TRG transcripts. Of note, CD8+ cytotoxic T cells were markedly less abundant than other αβ subsets (bottom) and were almost exclusively found in a single cat (22-007). (B) TRAV (top) and TRBV (bottom) gene usage of αβ T cell subsets. CD8+ cytotoxic T cells show differential gene usage compared to other αβ T cell subsets. The most abundant V genes are highlighted in either red (CD8+ cytotoxic) or blue (other αβ T cell subsets). (C) Junctional length and V/J usage of ab T cell subsets stratified by TRA/TRB vs. TRA/TRB/TRG expression. CD8+ cytotoxic T cells with TRA/TRB/TRG transcripts use distinct TRA and TRB but not TRG V/J gene combinations. TRG transcripts in this group are characterized by a limited junctional diversity. The most abundant V/J gene combinations in CD8+ cytotoxic T cells have been highlighted. (D) Position weight matrix of the most common TRA (top), TRB (middle) and TRG (bottom) junctional motifs in CD8+ cytotoxic T cells with TRA/TRB/TRG expression. (E) Combinatorial diversity of TRA and TRB chains in αβ T cells with TRA/TRB/TRG expression. Each stratum in the left and right axes represent a unique TRA V/J and TRB V/J combination, respectively. The connections between the left (TRA) and right (TRB) strata represent the pairing of specific TRA and TRB combinations. Compared to other subsets, CD8+ cytotoxic T cells (blue lines) have more frequent pairings between specific TRA and TRB combinations, suggesting that CD8+ cytotoxic T cells exhibit more focused TRA/TRB pairing patterns and lower combinatorial diversity. The three dominant pairings are: (1) TRAV23/TRAJ25 TRBV25/TRBJ2-6 (51 cells), (2) TRAV25/TRAJ41 TRBV25/TRBJ2-6 (35 cells), and (3) TRAV8-6/TRAJ30 TRBV4-2/TRBJ2-6 (27 cells).
The majority of γδ T cells expressed TRG/TRD transcripts (49.7%) followed by cells with an additional TRB transcripts (16.4%) (Figure 5A). Of note, 11.5% of cells classified as γδ cells based on their global transcriptome contained TRA/TRB/TRG transcripts, which might represent αβ T cells with an effector phenotype similar to γδ T cells. Next, we analyzed γδ T cells as defined by TCR expression rather than global gene expression. Out of 154 cells with a single productive TRD and TRG rearrangement each, 127 (82.5%) had been classified a γδ T cells based on their global transcriptome while 20 (13.0%) had been classified as CD4+ TEM cells. When assessing the pairing of TRD and TRG chains based on V gene usage, TRD V genes paired with TRG V genes in an all-vs-all pattern (Figure 6A). The exception to this rule were TRD chains utilizing the TRDV4 gene (n=19), which exclusively paired with chains containing the TRGV5-3 gene. TRDV4 and TRGV5-3 genes each rearranged to a single J gene only (TRDV4/TRDJ3 and TRGV5-3/TRGJ5-1). Interestingly, the V gene usage of cells coincided with functional phenotype as 18/19 of cells with TRDV4/TRGV5-3 usage were classified as CD4+ TEM based on their global transcriptome and only one was identified as γδ T cell. In addition, TRDV4/TRGV5-3 cells had a highly limited junctional length and sequence composition (Figures 6B, C). All cells had the same 17 aa TRD sequence ‘CASDIGGSSWDTRQMFF’ and 16/19 cells shared the 13 aa TRG sequence ‘CACWDESGWIKIF’. These cells were found in all four cats, suggesting that this might be a widely shared γδ T cell subset in cats.
Figure 6. Characterization of γδ T cells with one TRG & TRD rearrangement each. (A) Correlation of TRG/TRD V gene pairing and functional phenotype. Three TRD V genes (TRDV3, TRDV5-1, TRDV5-2, TRDV5-3) pair with three TRG V genes (TRGV2-1, TRGV2-2, TRGV2-4) in an all-vs-all fashion (white strata) while TRDV4 genes exclusively pair with TRDV5-3 genes (grey strata). Cells with TRDV4/TRDV5-3 gene usage have a CD4+ TEM phenotype (orange alluvium), while other cells have a predominant γδ phenotype (green alluvium). (B) Junctional diversity of γδ T cell subsets. TRDV4/TRDV5-3 γδ T cells have a highly skewed TRD & TRG junctional length and are found in all 4 cats. (C) Junctional sequence diversity of the TRD (top) & TRG (bottom) rearrangements. All TRDV4/TRDJ3 rearrangements and 16/19 TRG5-3/TRGJ5-1 rearrangements share the same junctional sequence, respectively. The three divergent TRG sequences were highly similar and either had a single amino acid variation in position 6 (2/3) or contained one additional histidine between positions 5 and 6 (CACWDHESGWIKIF, not shown).
Next, we characterized the repertoire overlap between cats and the relatedness of junctional regions across cats and T cell subsets (Figure 7). The number of shared clonotypes was highest for TRA and TRG but overall low (Figure 7A). When expanding the definition of shared repertoires to include clonotypes with one amino acid difference, the TRA locus had the highest number, size and publicity of clusters (Figures 7B, C).
Figure 7. (A) Shared clonotypes across 4 cats. The TRA locus exhibits the highest degree of publicity. The TRD and TRG clonotypes that are shared by all four cats are ‘CASDIGGSSWDTRQMFF’ and ‘CACWDESGWIKIF‘, respectively (see also Figure 6C) (B) Characterization of clusters based on mean centrality and cluster density. Higher values reflect larger and more connected clusters. (C) Network plots of TRA clusters with high centrality and density. Each node represents a unique clonotype, all clonotypes have identical junctional lengths and edges connect clonotypes with one amino acid sequence difference in the junctional region. The clusters contain clonotypes from all four cats and different T cell subsets. (D) Position weight matrices of TRA junctional regions of two representative clusters with high centrality and density (identified in Figure 7B).
2.6 Feline circulating monocytic cells
Although CD5+ enrichment was performed for T cell selectivity, we serendipitously captured neutrophils and monocytic cells, which provides insight into the transcriptome of circulating myeloid cells in cats. Clustering analysis of both cell types was performed independently. Due to the low myeloid cell count in the 1YO, it was excluded from this analysis.
Monocytic cells clustered into 5 groups and were largely represented by the 9YO sample (Figure 8A; Supplementary Presentation 1). Cluster 3 represented neutrophils which is likely due to the close clustering of myeloid cells on the global object (Figure 1B; Supplementary Presentation 1). Cluster 0, 1 and 2 were identified as classical monocytes (CM) with all clusters having high expression of markers PLBD1, CD83 and classical monocyte markers CD14 and VCAN (28, 29, 46) (Figures 8B–E). We did not identify any non-classical monocytes, based on absence of CD16 expression, as seen in other species such as humans, horses and cattle (22, 29, 46, 47). Cluster 4 showed an upregulation of receptor-type tyrosine-protein kinase FLT3, a conventional pan-dendritic cell marker (48) (Figure 8F). In addition, these cells exhibited upregulation of MHC class-II molecules (FECA-DRB, HLA-DRA, HLA-DOA, and HLA-DMD) and were hence identified as conventional dendritic cells (cDC) (49) (Figure 8G).
Figure 8. Feline circulating monocytes include classical monocytes with more differentiated clusters with traditional proinflammatory and unconventional anti-inflammatory phenotypes. Unsupervised clustering revealed 5 clusters. (A) UMAP of scRNA-seq atlas of monocytes. (B–F) UMAP of monocytes colored by expression levels of canonical markers. (G) Dot plot of genes of monocytic phenotypes across the 4 monocytic clusters. (H–J) Dot plots of top GO biological process terms called based on sets of positive differentially expressed genes identified via Seurat FindMarker function (Wilcoxon rank sum, Adj P <0.05). CM, Classical monocyte; cDC, Conventional dendritic cell; Pro-inf, proinflammatory; Anti-inf, anti-inflammatory..
In-depth assessment of differentially expressed genes across the four monocytic clusters identified signatures that further characterized these cell types. Cluster 2 demonstrated a proinflammatory phenotype characterized by upregulation of classic proinflammatory cytokines (IL1A, IL1B), TNF/NF-KB pathway (TNF, TNFAIP3, NFKBIA, NFKBIZ, NFKB1) and M1-macrophage-associated proinflammatory cytokines (CXCL10, CXCL16, CCL3, CCL4 and CCL5) (Figures 8G, H; Supplementary Tables 6, 7) (20, 50–54). Cluster 1 gene ontology revealed an upregulation of leukocyte migration and chemotaxis terms in addition to upregulation of proinflammatory cytokines IL1A and IL1B (Figures 8G, I; Supplementary Tables 6, 7). However, these cells also lack TNF expression and upregulate anti-inflammatory cytokines IL10 and IL18 and proteins which promote anti-inflammatory phenotypes in macrophages such as LCN2 and MSRB1 (55, 56) (Figure 8G). They have the lowest expression of CD86, a T cell co-stimulatory molecule of antigen presenting cells, which is associated with an anti-inflammatory phenotype (57). Most strikingly, these cells were found to upregulate genes similar to a novel monocyte cluster found in SARS-CoV-2 infected and recovered patients. These genes include AREG, EREG, BCL6, IL10, and IL18 which have been described as anti-inflammatory tissue repair genes (Figure 8G) (58, 59). Thus, these cells were designated as anti-inflammatory, although they do not follow the traditional M2 polarization phenotype described in human monocytes (60). Our data suggests a novel classical monocyte subtype with an anti-inflammatory phenotype.
Cluster 0 (CM) appears the least differentiated with low expression of inflammatory cytokines and the top GO terms include more cellular rather than immune related functions (Figure 8H; Supplementary Table 7). Cluster 4 (cDC) showed an upregulation of and the GO terms associated with leukocyte activation and T cell activation via TCR contact with antigen bound MHC (Figures 8G, J; Supplementary Tables 6, 7). Given the strong antigen presentation phenotype, the annotation of cDC was further corroborated.
2.7 Feline circulating neutrophils
We identified 2 major neutrophil clusters (Figure 9A) based on the expression of canonical neutrophil markers including CSF3R and ELANE as well as mature neutrophil functional markers vimentin (VIM), thioredoxin (TXN), S100A9 and S100A12, suggesting both clusters are comprised of mature neutrophils (25, 61, 62) (Figures 9B–G). Cluster 1 exhibited several notable differences compared to cluster 0. Differential gene expression analysis accompanied by gene ontology showed an upregulation in cluster 1 of apoptotic signaling and activated neutrophil functions such as response to lipopolysaccharide and IL-2 production (Figure 9H; Supplementary Tables 8, 9). Thus, we annotated this cluster as activated neutrophils. When investigating top differentially expressed genes, we noted a higher expression of IL1A and TNFAIP3, suggesting a proinflammatory phenotype (63) (Figures 9I, J). Similarly, cluster 1 had an increased expression of VCAN, which is upregulated in the skin in response to damage from ultraviolet light B and reactive oxygen species (64) (Figure 9K). Lastly, cluster 1 exhibited lower RNA counts per cell, which is consistent with findings in murine neutrophils where the number of genes and counts of RNA decrease with neutrophil differentiation, maturation and eventual activation and function (65) (Figures 9L, M). Since a pathway driving this clustering could not be resolved, we examined an interferon (IFN) module score, which is a composite score based on expression levels of interferon stimulated genes (Figure 9N). Cluster 1 did not show an increased IFN score but the adjacent cells in cluster 0 did (Figure 9N). Interferon stimulation can delay neutrophil apoptosis, which could represent one of the potential driving factors for cells to move from cluster 0 to cluster 1 (66). This analysis is the first to resolve circulating neutrophil transcriptomic signatures in cats at the single cell level.
Figure 9. Feline circulating neutrophils separate into 2 clusters based on activation state. (A) UMAP of scRNA-seq atlas of neutrophils. (B–G) UMAP of neutrophils colored by expression levels of canonical and functional neutrophil markers. (H) Dot plots of top GO biological process terms in cluster 1 called based on sets of positive differentially expressed genes identified via Seurat FindMarker function (Wilcoxon rank sum, Adj P <0.05). (I–K) UMAP of neutrophils colored by expression levels of select differentially upregulated genes in cluster 1. (L, M) Violin plots of RNA counts and feature (unique RNA) counts per cell by cluster; (****) indicates P value less than 0.001 by T-test. (N) UMAP of neutrophils colored by interferon (IFN) gene composite score per cell; genes included are named in figure.
3 Discussion
The importance of the cat as a companion animal as well as a model for infectious diseases warrants an in-depth characterization and understanding of the feline immune system. Our study is one of the first to characterize feline immune cell subpopulations at the single cell level. We utilized 5’ single cell RNA-sequencing with V(D)J analysis to resolve the heterogeneity of CD5+ enriched peripheral blood immune cells. We resolved populations of T cells, neutrophils, B cells and monocytes. Additionally, this atlas is the first to annotate the single cell transcriptome of rarer cell types, plasmacytoid dendritic cells and mast cells, which are seldomly captured in single cell atlases (67, 68).
Among T cells, we were able to resolve different populations of naive cells including a significant number of CD8A/CD8B naive cells. We further identified T effector cells and terminally differentiated cytotoxic effectors. The transcriptomic gradient of T cell differentiation was found to fit the paradigm seen with other species (23). Additionally, we captured feline γδ T cells for the first time. Upon investigation of effector cells, we were able to demonstrate known T helper effector phenotypes among the most differentiated effectors. Of note, the Treg phenotype was found to be the most prominent. The importance of Treg populations has been studied in cats, particularly in the context of feline mammary carcinoma (69–71). Our study corroborates the large presence of this phenotype and the potential importance of Treg cells in feline immunology.
We performed a cross species integrative analysis of T cells to take advantage of the growing database of single cell RNA data. We found relatively high conservation of T cell subtypes along an effector gradient with relatively equitable representation of other veterinary species (horse, dog, pig) and humans with the cat. This analysis also revealed the potential negative impacts of CD5+ enrichment for T cells. Effector clusters with high CD5 expression were found to be underrepresented in the cat. Thus, a broader T cell selection protocol would be necessary in future studies to further investigate effector cells in cats.
In addition to characterizing T cells based on their transcriptome, we analyzed TCR expression on a single cell level for the first time in the cat. Cytotoxic CD8+ T cells differed from other αβ T cell subsets in multiple aspects. In addition to the expected TRA and TRB transcripts, a high percentage of cytotoxic CD8+ T cells expressed a productive TRG transcript. Given that this was not observed in naïve CD8+ T cells, we attribute this phenomenon to differential transcriptional regulation of naive and effector T cells. Compared to other αβ T cell subsets, cytotoxic T cells also exhibited divergent TRAV and TRBV gene usage but not TRGV gene usage. Moreover, the junctional length of TRG transcripts was skewed towards a highly conserved 16 aa motif. This finding is of potential significance for clonality testing, an adjunct method for differentiating reactive from neoplastic lymphoid proliferations. If the majority of cytotoxic CD8+ T cells harbor a length restricted TRG rearrangement, then reactive CD8+ dominant T cell proliferations might exhibit reduced TRG diversity, which could result in a higher false positive rate. Of note, since other αβ T cell subsets did not express a significant number of productive TRG transcripts, it is unknown if similar length restricted rearrangements occur in these subsets on a DNA level. A previous study examining the feline TCR transcriptome across thymus, spleen and mesenteric lymph nodes did not identify any TRG length restrictions (13). However, the previous study aimed to identify global rearrangement patterns and hence utilized a bulk sequencing strategy, which is not suited to resolve repertoire characteristics of individual T cell subsets.
Cytotoxic CD8+ T cells with TRG transcripts utilized few but distinct TRA and TRB V/J parings and had conserved junctional motifs. A previous study utilizing scRNA-seq with VDJ sequencing in dogs identified MAIT-like cells with restricted TRAV gene usage and junctional length as well as numerous distinct but highly similar TRA gene sequences (19). A cluster analysis of TCR sequences in this study did not identify T cells with similarly characteristics. The fact that we were unable to identify MAIT-like cells in cats could reflect their true absence or represent an artifact of CD5-based enrichment, which might have resulted in the exclusion of CD5- or CD5low subsets. Additional studies with more animals, inclusion of different tissue types and without selection bias are needed to elucidate the existence of innate-like αβ T cells in cats.
The annotation of γδ T cells based on their TCR expression rather than global transcriptome revealed a population of γδ T cells with invariant TRD and TRG V/J gene usage and a highly conserved junctional region. The existence of T cells with TRGV5-3 gene usage and reduced TCR diversity had previously been suspected based on the finding of a recurrent clone in liver and intestinal samples of cats (87). Our results demonstrate that this canonical rearrangement stems from γδ T cells with a similarly restricted TRD chain utilizing TRDV4 and a highly conserved junctional motif. This population was missed during an initial analysis because cells clustered with and were obscured by more numerous CD4+ αβ TEM cells when annotated by global transcriptome. Akin to CD8+ αβ T cells with a canonical TRG rearrangement, TRGV5-3/TRDV4 γδ T cells could result in false positive clonality testing results and hence require careful primer set design and result interpretation.
T cell lineage annotation based on global transcriptome did not always match the classification based on T cell receptor expression. In addition to the previously discussed TRGV5-3/TRDV4 γδ T cells that masqueraded as CD4+ αβ TEM cells, annotation based on global transcriptome revealed presumed γδ T cells with TRA-TRB and TRA-TRB-TRG rearrangements. Given that TCR cell surface expression is the ultimate arbiter for T cell lineage, these cells most likely represented αβ T cells with functional properties similar to that of γδ T cells. Additional factors that might result in erroneous lineage classification are low TCR expression or incomplete annotation of TCR genes in the cat reference transcriptome. In mouse models it has been demonstrated that adaptive γδ T cells that have overlapping phenotypes with the αβ T cells and it has been suggested that CD8+ γδ T cells may be functionally indistinguishable from CD8+ αβ T cells (88). Consequently, αβ vs. γδ lineage classification based on global transcriptome alone might not always be always accurate and additional studies with concomitant analysis of all four TCR loci are needed to further substantiate this phenomenon.
Among the myeloid cells, we resolved 3 clusters of classical monocytes and 1 cluster of conventional dendritic cells. Of the three clusters of monocytes, we found two to be differentiated in a somewhat polarized manner- pro-inflammatory versus anti-inflammatory. The anti-inflammatory population was unconventional in that it had an upregulation of classic pro-inflammatory cytokines IL-1A and IL-1B but also an upregulation of genes associated with a unique and novel monocytic cluster found in human COVID-19 patients such as AREG, EREG, IL10, IL18 and BCL6 (59, 72). Sub clustering of neutrophils identified a larger mature cluster and a smaller exhausted/activated cluster which showed upregulation of pro-inflammatory and apoptotic genes. Thus, analysis of these myeloid cells revealed the presence of conventional and non-conventional subtypes, providing a footing for feline circulating myeloid immunology.
Currently, feline T cell subtypes have only been described at the flow cytometric level, and with a focus on CD4/CD8 populations due to paucity of species-specific reagents. Our data indicates the presence of unique subtypes of peripheral blood immune cells including αβ CD8+ cytotoxic cells, CD8A/CD8B naive T cells, plasmacytoid DC, mast cells, conventional DCs, anti-inflammatory monocytes and mature and end stage neutrophils. Our study is the first of many necessary for the further characterization of the feline immune system. Taken together, it has increased our understanding of feline circulating immune cells, including T cell subtyping and capture of rare populations previously undescribed in the cat.
4 Materials and methods
4.1 Samples and Sorting
Whole blood was obtained from four healthy, female, domestic shorthair cats (13-101/9 years, 18-003/6 years, 21-005/1 year, 22-007/6 months) housed at the UC Davis Nutrition and Pet Care Center in accordance with the UC Davis Policy and Procedure Manual section 290-30 and Animal Use and Care protocol #22150. Blood was collected in purple top EDTA tubes and immediately processed in one batch without freezing until completion of the reverse transcription step during library preparation. Peripheral blood mononuclear cells (PBMC) were isolated from whole blood by density gradient centrifugation. Briefly, 5 mL of whole blood was layered on top of 10 mL Histopaque 1077 (density 1.077g/mL, Sigma-Aldrich, Burlington, MA, USA) in a 50 ml tube. The samples were centrifuged at 145 relative centrifugal force (RCF) for 30 minutes at room temperature (RT) without break. After centrifugation, the PBMC layer at the plasma-histopaque interface was aspirated and transferred to a clean 15 ml tube. The PBMCs were then washed twice with 14 ml PBMC isolation media (500 mL Hanks’ Balanced Salt Solution (HBSS), 15 mL Fetal Bovine Serum (FBS), 2 mL EDTA) by centrifuging at 850 RCF for 10 minutes at RT. Following the second washing step, the PBMC pellet was resuspended in 200 μl modified flow buffer solution (MFB, 500 ml phosphate buffer saline (PBS), 0.467g EDTA tetrasodium salt and 1% Horse serum).
The isolated PBMCs were first incubated with a primary anti-feline CD5 monoclonal antibody (clone FE 1.1.B11-IgG1-Azide, Peter F. Moore Leukocyte Antigen Biology Laboratory, UC Davis, CA, USA) or a non-cross reacting canine antibody (negative control) for 30 min at RT. This was followed by a wash step (addition off 3 mL of MFB, followed by centrifugation at 400 RCF for 3 min at RT and resuspension in 100 - 200 μl of supernatant), and incubation with a secondary horse-anti-mouse FITC antibody (1:99 in MFB) for 15 minutes at RT in darkness. Finally, the stained PBMC was washed again (as above), resuspended in 500 μl of MFB and stained with DAPI (1 μg/mL). Cells were sorted at the UC Davis Comprehensive Cancer Center Flow Cytometry Core. Cells were initially gated based on FSC vs. SSC to exclude debris and aggregates. Singlet gating was then applied based on SSC to exclude doublets, followed by discrimination of live and dead cells based on DAPI staining and subsequent gating for CD5 positive (FITC positive) cells (Supplementary Presentation 2).
4.2 Single-cell 5’ RNA and V(D)J library preparation, sequencing and mapping
Libraries were constructed using the Chromium Next GEM Single Cell 5’ Kit v2 (10x Genomics). All four feline T cell receptor loci were amplified using the Chromium Single Cell V(D)J Enrichment Kit according to the manufacturer’s instructions but with custom reverse primers specific for the feline orthologs (Supplementary Table 10). Separate libraries were prepared for the transcriptome, αβ TCRs, and γδ TCRs (Supplementary Table 10). Libraries were sequenced using the Illumina NovaSeq S4 platform using 150 paired-end reads to a target read depth of 30,000 reads per cell. Raw FASTQ sequencing data were processed by the UC Davis Bioinformatics Core using 10X CellRanger version 7.1 for alignment to the cat reference genome (Felis_catus_9.0). To improve the quality of annotation of the Felis catus genome, we converted unmapped ENSEMBL IDs to human homologs from BioMart for greater interpretative power. Gene mapping was performed in the following order: (1) Mapping of one-to-one orthologs from cat to human (2) Mapping of one-to-many and many-to-many was performed by choosing a representative cat gene based on the following parameters, in respective order [orthology_confidence, gene_order_confidence, Whole_Genome_alignment, %_query_identical_target, %target_identical_query]. Additionally, specific unmapped genes were manually annotated or filtered (Supplementary Table 11). Further analysis was then continued using Seurat v5 (73).
4.3 Count matrix pre-processing and quality control
Filtered feature matrices for the 4 samples were independently normalized and scaled using the function SCTransform v2 (glmGamPoi) (74). Dimensionality reduction (determined using knee identification in elbow plot), k-means clustering (resolution determined using clustree package) and neighbor identification (74–76) was performed. Low quality clusters were identified and filtered based on low gene counts, high mitochondrial read percentage and lack of biologically relevant markers. To further ensure quality, doublet identification was performed using package DoubletFinder according to package vignette and ambient RNA decontamination was performed using the DecontX package as per package tutorials adapted to the data (77, 78).
4.4 Integration and downstream analysis
Objects were integrated by sample using the Seurat native reciprocal principal component analysis (RPCA) (79). The embeddings were utilized for Uniform Manifold Approximation and Projection (UMAP) and subsequent clustering and neighbor identification as described prior. Differential gene expression analysis was performed using each cluster identified by the Seurat function FindAllMarkers (Non-parametric Wilcoxon Rank Sum with adjusted P value <0.05). Gene set enrichment analyses including over-enrichment analysis (ORA) of gene ontology (biological process) was performed using ClusterProfiler on human homologized genes (adjusted P value <0.05) (80). Additional visualizations and associated statistics were performed using ggplot2 and gg supplementary packages (81). Figures were illustrated on Bio Render.
4.5 Cell type sub clustering
Cell type sub-clustering (example: T cell sub-clustering) analysis was performed in the same manner as described prior. Pseudotime values were calculated using the UMAP embeddings produced via Seurat with RPCA integration using package monocle3 (82). Gene module scores were calculated via the Seurat function AddModuleScore.
4.6 Cross species T cell integrative analysis
Published single cell RNA-sequencing datasets were obtained from NCBI GEO Accession for 4 species including dog [Canis familiaris: GSE225599 (18)], human (Homo sapiens: https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/pbmc_10k_v3), horse [Equus caballus: GSE148416 (22)] and pig [Sus scrofa: https://data.nal.usda.gov/dataset/data-reference-transcriptomics-porcine-peripheral-immune-cells-created-through-bulk-and-single-cell-rna-sequencing (21)]. All species genes were mapped to human homologs as described prior for the cat. Only homologized shared genes between all 5 species were utilized for analysis. T cell clusters annotated by authors were validated using canonical T cell markers (CD3E, CD3D, TCF7, LCK, ITK). T cells were integrated across species using Seurat canonical correlation analysis (CCA) and downstream analysis proceeded in the same manner as described prior for global object analysis.
4.7 VDJ analysis
Contigs were mapped to available IMGT reference genes (as of February 2023) using Cell Ranger version 7.1. All analyses were based on the ‘all_contig_annotations.csv’ file using the R packages tidyverse, circlize and ggseqlogo (83–85). A clonotype was defined as a unique TCR amino acid sequence. To exclude incomplete contigs or inter-locus chimeras, we only retained productive contigs for which V, J and C genes could be identified and belonged to the same locus. Cluster analysis was performed as previously described (86).
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/, GSE267355.
Ethics statement
The animal study was approved by UC Davis' Institutional Animal Care and Use Committee in accordance with the UC Davis Policy and Procedure Manual section 290-30 and Animal Use and Care protocol #22150. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
RR: Validation, Investigation, Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Formal analysis, Data curation. JW: Investigation, Writing – review & editing, Writing – original draft, Methodology, Funding acquisition, Formal analysis. HC: Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Data curation. PM: Writing – review & editing, Writing – original draft, Resources, Methodology, Investigation, Funding acquisition, Conceptualization. WV: Writing – review & editing, Writing – original draft, Methodology, Investigation, Funding acquisition, Conceptualization. SK: Visualization, Validation, Supervision, Software, Resources, Project administration, Formal analysis, Data curation, Writing – review & editing, Writing – original draft, Methodology, Investigation, Funding acquisition, Conceptualization.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This project was funded by grant #2022-24-F from the Center for Companion Animal Health, University of California, Davis. The project was supported by the University of California Davis Flow Cytometry Shared Resource Laboratory with funding from the NCI P30 CA093373 (Comprehensive Cancer Center) and S10 OD018223, with technical assistance from Bridget McLaughlin, Jonathan Van Dyke and Ashley Karajeh. The sequencing was carried out by the DNA Technologies and Expression Analysis Core at the UC Davis Genome Center, supported by NIH Shared Instrumentation Grant 1S10OD010786-01. Publication fees were funded by the UC Davis Library Open Access Fund and the Center for Companion Animal Health, UC Davis.
Acknowledgments
We thank Jie (Jessie) Li from the UC Davis Bioinformatics Core for data analysis and Kristy Harmon from the Leukocyte Biology Antigen Laboratory for her flow cytometry support.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1438004/full#supplementary-material
References
1. Rowan AN. Companion Animal Statistics in the USA. (2018). Available online at: https://www.wellbeingintlstudiesrepository.org/demscapop/7/. (Accessed February 20, 2024).
2. Larkin M, Radich R, AVMA. Am Vet Med Assoc. Pet population still on the rise, with fewer pets per household (2021). Available online at: https://avmajournals.avma.org/display/post/news/pet-population-still-on-the-rise–with-fewer-pets-per-household.xml?print&print&print&print&print. (Accessed February 20, 2024).
3. Fromont E, Artois M, Langlais M, Courchamp F, Pontier D. Modelling the feline leukemia virus (FeLV) in natural populations of cats (Felis catus). Theor Popul Biol. (1997) 52:60–70. doi: 10.1006/tpbi.1997.1320
4. Powers JA, Chiu ES, Kraberger SJ, Roelke-Parker M, Lowery I, Erbeck K, et al. Feline leukemia virus (FeLV) disease outcomes in a domestic cat breeding colony: relationship to endogenous feLV and other chronic viral infections. J Virol. (2018) 92. doi: 10.1128/JVI.00649-18
5. Helfer-Hungerbuehler AK, Shah J, Meili T, Boenzli E, Li P, Hofmann-Lehmann R. Adeno-associated vector-delivered CRISPR/saCas9 system reduces feline leukemia virus production in vitro. Viruses. (2021) 13. doi: 10.3390/v13081636
6. Bienzle D. FIV in cats–a useful model of HIV in people? Vet Immunol Immunopathol. (2014) 159:171–9. doi: 10.1016/j.vetimm.2014.02.014
7. Paillot R, Richard S, Bloas F, Piras F, Poulet H, Brunet S, et al. Toward a detailed characterization of feline immunodeficiency virus-specific T cell immune responses and mediated immune disorders. Vet Immunol Immunopathol. (2005) 106:1–14. doi: 10.1016/j.vetimm.2004.12.023
8. Endo Y, Matsumura S, Washizu T, Ishida T. Alteration of T-cell subsets in the lymph nodes from cats infected with feline immunodeficiency virus. J Vet Med Sci. (1997) 59:739–46. doi: 10.1292/jvms.59.739
9. Hoffmann-Fezer G, Mortelbauer W, Hartmann K, Mysliwietz J, Thefeld S, Beer B, et al. Comparison of T-cell subpopulations in cats naturally infected with feline leukaemia virus or feline immunodeficiency virus. Res Vet Sci. (1996) 61:222–6. doi: 10.1016/s0034-5288(96)90067-3
10. Day MJ. Cats are not small dogs: is there an immunological explanation for why cats are less affected by arthropod-borne disease than dogs? Parasit Vectors. (2016) 9:507. doi: 10.1186/s13071-016-1798-5
11. Akira S, Uematsu S, Takeuchi O. Pathogen recognition and innate immunity. Cell. (2006) 124:783–801. doi: 10.1016/j.cell.2006.02.015
12. Jiao A, Zhang C, Wang X, Sun L, Liu H, Su Y, et al. Single-cell sequencing reveals the evolution of immune molecules across multiple vertebrate species. J Advert Res. (2024) 55:73–87. doi: 10.1016/j.jare.2023.02.017
13. Radtanakatikanon A, Keller SM, Darzentas N, Moore PF, Folch G, Nguefack Ngoune V, et al. Topology and expressed repertoire of the Felis catus T cell receptor loci. BMC Genomics. (2020) 21:20. doi: 10.1186/s12864-019-6431-5
14. Ginhoux F, Yalin A, Dutertre CA, Amit I. Single-cell immunology: Past, present, and future. Immunity. (2022) 55:393–404. doi: 10.1016/j.immuni.2022.02.006
15. Jovic D, Liang X, Zeng H, Lin L, Xu F, Luo Y. Single-cell RNA sequencing technologies and applications: A brief overview. Clin Transl Med. (2022) 12:e694. doi: 10.1002/ctm2.694
16. Sun X, Lin X, Li Z, Wu H. A comprehensive comparison of supervised and unsupervised methods for cell type identification in single-cell RNA-seq. Brief Bioinform. (2022) 23. doi: 10.1093/bib/bbab567
17. Wang M, Song WM, Ming C, Wang Q, Zhou X, Xu P, et al. Guidelines for bioinformatics of single-cell sequencing data analysis in Alzheimer’s disease: review, recommendation, implementation and application. Mol Neurodegener. (2022) 17:17. doi: 10.1186/s13024-022-00517-z
18. Ammons DT, Harris RA, Hopkins LS, Kurihara J, Weishaar K, Dow S. A single-cell RNA sequencing atlas of circulating leukocytes from healthy and osteosarcoma affected dogs. Front Immunol. (2023) 14:1162700. doi: 10.3389/fimmu.2023.1162700
19. Eschke M, Moore PF, Chang H, Alber G, Keller SM. Canine peripheral blood TCRαβ T cell atlas: Identification of diverse subsets including CD8A+ MAIT-like cells by combined single-cell transcriptome and V(D)J repertoire analysis. Front Immunol. (2023) 14:1123366. doi: 10.3389/fimmu.2023.1123366
20. Gao Y, Li J, Cai G, Wang Y, Yang W, Li Y, et al. Single-cell transcriptomic and chromatin accessibility analyses of dairy cattle peripheral blood mononuclear cells and their responses to lipopolysaccharide. BMC Genomics. (2022) 23:338. doi: 10.1186/s12864-022-08562-0
21. Herrera-Uribe J, Wiarda JE, Sivasankaran SK, Daharsh L, Liu H, Byrne KA, et al. Reference transcriptomes of porcine peripheral immune cells created through bulk and single-cell RNA sequencing. Front Genet. (2021) 12:689406. doi: 10.3389/fgene.2021.689406
22. Patel RS, Tomlinson JE, Divers TJ, Van de Walle GR, Rosenberg BR. Single-cell resolution landscape of equine peripheral blood mononuclear cells reveals diverse cell types including T-bet+ B cells. BMC Biol. (2021) 19:13. doi: 10.1186/s12915-020-00947-5
23. Li Z, Sun C, Wang F, Wang X, Zhu J, Luo L, et al. Molecular mechanisms governing circulating immune cell heterogeneity across different species revealed by single-cell sequencing. Clin Transl Med. (2022) 12:e689. doi: 10.1002/ctm2.689
24. Touw IP, Palande K, Beekman R. Granulocyte colony-stimulating factor receptor signaling: implications for G-CSF responses and leukemic progression in severe congenital neutropenia. Hematol Oncol Clin North Am. (2013) 27:61–73, viii. doi: 10.1016/j.hoc.2012.10.002
25. Sprenkeler EGG, Zandstra J, van Kleef ND, Goetschalckx I, Verstegen B, Aarts CEM, et al. S100A8/A9 is a marker for the release of neutrophil extracellular traps and induces neutrophil activation. Cells. (2022) 11. doi: 10.3390/cells11020236
26. Akgun E, Tuzuner MB, Sahin B, Kilercik M, Kulah C, Cakiroglu HN, et al. Proteins associated with neutrophil degranulation are upregulated in nasopharyngeal swabs from SARS-CoV-2 patients. PloS One. (2020) 15:e0240012. doi: 10.1371/journal.pone.0240012
27. Blanter M, Cambier S, De Bondt M, Vanbrabant L, Pörtner N, Abouelasrar Salama S, et al. Method matters: effect of purification technology on neutrophil phenotype and function. Front Immunol. (2022) 13:820058. doi: 10.3389/fimmu.2022.820058
28. Buechler C, Ritter M, Orsó E, Langmann T, Klucken J, Schmitz G. Regulation of scavenger receptor CD163 expression in human monocytes and macrophages by pro- and antiinflammatory stimuli. J Leukoc Biol. (2000) 67:97–103. https://www.ncbi.nlm.nih.gov/pubmed/10648003.
29. Sampath P, Moideen K, Ranganathan UD, Bethunaickan R. Monocyte subsets: phenotypes and function in tuberculosis infection. Front Immunol. (2018) 9:1726. doi: 10.3389/fimmu.2018.01726
30. Reshetnikova E, Guselnikov S, Volkova O, Baranov K, Taranin A, Mechetina L. B cell-specific protein FCRLA is expressed by plasmacytoid dendritic cells in humans. Cytometry B Clin Cytom. (2018) 94:683–7. doi: 10.1002/cyto.b.21611
31. Rissoan MC, Duhen T, Bridon JM, Bendriss-Vermare N, Péronne C, de Saint Vis B, et al. Subtractive hybridization reveals the expression of immunoglobulin-like transcript 7, Eph-B1, granzyme B, and 3 novel transcripts in human plasmacytoid dendritic cells. Blood. (2002) 100:3295–303. doi: 10.1182/blood-2002-02-0638
32. Edwards AD, Diebold SS, Slack EMC, Tomizawa H, Hemmi H, Kaisho T, et al. Toll-like receptor expression in murine DC subsets: lack of TLR7 expression by CD8α+ DC correlates with unresponsiveness to imidazoquinolines. Eur J Immunol. (2003) 33:827–33. doi: 10.1002/eji.200323797
33. Medina KL, Tangen SN, Seaburg LM, Thapa P, Gwin KA, Shapiro VS. Separation of plasmacytoid dendritic cells from B-cell-biased lymphoid progenitor (BLP) and Pre-pro B cells using PDCA-1. PloS One. (2013) 8:e78408. doi: 10.1371/journal.pone.0078408
34. Arroyo Hornero R, Idoyaga J. Plasmacytoid dendritic cells: A dendritic cell in disguise. Mol Immunol. (2023) 159:38–45. doi: 10.1016/j.molimm.2023.05.007
35. Atiakshin D, Kostin A, Trotsenko I, Samoilova V, Buchwalow I, Tiemann M. Carboxypeptidase A3-A key component of the protease phenotype of mast cells. Cells. (2022) 11(3):570. doi: 10.3390/cells11030570
36. Cruse G, Kaur D, Leyland M, Bradding P. A novel FcϵRIβ-chain truncation regulates human mast cell proliferation and survival. FASEB J. (2010) 24:4047–57. doi: 10.1096/fj.10-158378
37. Lee H, Joo JY, Kang J, Yu Y, Kim YH, Park HR. Single-cell analysis of platelets from patients with periodontitis and diabetes. Res Pract Thromb Haemost. (2023) 7:100099. doi: 10.1016/j.rpth.2023.100099
38. Yang S, Liu F, Wang QJ, Rosenberg SA, Morgan RA. The shedding of CD62L (L-selectin) regulates the acquisition of lytic activity in human tumor reactive T lymphocytes. PloS One. (2011) 6:e22560. doi: 10.1371/journal.pone.0022560
39. Ivetic A, Hoskins Green HL, Hart SJ. L-selectin: A major regulator of leukocyte adhesion, migration and signaling. Front Immunol. (2019) 10:1068. doi: 10.3389/fimmu.2019.01068
40. Bjorkdahl O, Barber KA, Brett SJ, Daly MG, Plumpton C, Elshourbagy NA, et al. Characterization of CC-chemokine receptor 7 expression on murine T cells in lymphoid tissues. Immunology. (2003) 110:170–9. doi: 10.1046/j.1365-2567.2003.01727.x
41. Cano-Gamez E, Soskic B, Roumeliotis TI, So E, Smyth DJ, Baldrighi M, et al. Single-cell transcriptomics identifies an effectorness gradient shaping the response of CD4+ T cells to cytokines. Nat Commun. (2020) 11:1801. doi: 10.1038/s41467-020-15543-y
42. Zhu J, Yamane H, Paul WE. Differentiation of effector CD4 T cell populations*. Annu Rev Immunol. (2009) 28:445–89. doi: 10.1146/annurev-immunol-030409-101212
43. Andreatta M, Corria-Osorio J, Müller S, Cubas R, Coukos G, Carmona SJ. Interpretation of T cell states from single-cell transcriptomics data using reference atlases. Nat Commun. (2021) 12:2965. doi: 10.1038/s41467-021-23324-4
44. Wang X, Shen X, Chen S, Liu H, Hong N, Zhong H, et al. Reinvestigation of classic T cell subsets and identification of novel cell subpopulations by single-cell RNA sequencing. J Immunol. (2022) 208:396–406. doi: 10.4049/jimmunol.2100581
45. Vallejo J, Saigusa R, Gulati R, Armstrong Suthahar SS, Suryawanshi V, Alimadadi A, et al. Combined protein and transcript single-cell RNA sequencing in human peripheral blood mononuclear cells. BMC Biol. (2022) 20:193. doi: 10.1186/s12915-022-01382-4
46. Schneider L, Marcondes NA, Hax V, da Silva Moreira IF, Ueda CY, Piovesan RR, et al. Flow cytometry evaluation of CD14/CD16 monocyte subpopulations in systemic sclerosis patients: a cross sectional controlled study. Adv Rheumatol. (2021) 61:27. doi: 10.1186/s42358-021-00182-8
47. Barut GT, Kreuzer M, Bruggmann R, Summerfield A, Talker SC. Single-cell transcriptomics reveals striking heterogeneity and functional organization of dendritic and monocytic cells in the bovine mesenteric lymph node. Front Immunol. (2022) 13:1099357. doi: 10.3389/fimmu.2022.1099357
48. Waskow C, Liu K, Darrasse-Jèze G, Guermonprez P, Ginhoux F, Merad M, et al. The receptor tyrosine kinase Flt3 is required for dendritic cell development in peripheral lymphoid tissues. Nat Immunol. (2008) 9:676–83. doi: 10.1038/ni.1615
49. ten Broeke T, Wubbolts R, Stoorvogel W. MHC class II antigen presentation by dendritic cells regulated through endosomal sorting. Cold Spring Harb Perspect Biol. (2013) 5:a016873. doi: 10.1101/cshperspect.a016873
50. Eislmayr K, Bestehorn A, Morelli L, Borroni M, Vande Walle L, Lamkanfi M, et al. Nonredundancy of IL-1α and IL-1β is defined by distinct regulation of tissues orchestrating resistance versus tolerance to infection. Sci Adv. (2022) 8:eabj7293. doi: 10.1126/sciadv.abj7293
51. Lopez-Castejon G, Brough D. Understanding the mechanism of IL-1β secretion. Cytokine Growth Factor Rev. (2011) 22:189–95. doi: 10.1016/j.cytogfr.2011.10.001
52. Hadadi E, Zhang B, Baidžajevas K, Yusof N, Puan KJ, Ong SM, et al. Differential IL-1β secretion by monocyte subsets is regulated by Hsp27 through modulating mRNA stability. Sci Rep. (2016) 6:39035. doi: 10.1038/srep39035
53. Wan J, Shan Y, Fan Y, Fan C, Chen S, Sun J, et al. NF-κB inhibition attenuates LPS-induced TLR4 activation in monocyte cells. Mol Med Rep. (2016) 14:4505–10. doi: 10.3892/mmr.2016.5825
54. Mussbacher M, Derler M, Basílio J, Schmid JA. NF-κB in monocytes and macrophages - an inflammatory master regulator in multitalented immune cells. Front Immunol. (2023) 14:1134661. doi: 10.3389/fimmu.2023.1134661
55. Guo H, Jin D, Chen X. Lipocalin 2 is a regulator of macrophage polarization and NF-κB/STAT3 pathway activation. Mol Endocrinol. (2014) 28:1616–28. doi: 10.1210/me.2014-1092
56. Lee BC, Lee SG, Choo MK, Kim JH, Lee HM, Kim S, et al. Selenoprotein MsrB1 promotes anti-inflammatory cytokine gene expression in macrophages and controls immune response in vivo. Sci Rep. (2017) 7:5119. doi: 10.1038/s41598-017-05230-2
57. Vishnyakova P, Poltavets A, Karpulevich E, Maznina A, Vtorushina V, Mikhaleva L, et al. The response of two polar monocyte subsets to inflammation. BioMed Pharmacother. (2021) 139:111614. doi: 10.1016/j.biopha.2021.111614
58. Zhang Y, Wang S, Xia H, Guo J, He K, Huang C, et al. Identification of monocytes associated with severe COVID-19 in the PBMCs of severely infected patients through single-cell transcriptome sequencing. Engineering. (2022) 17:161–9. doi: 10.1016/j.eng.2021.05.009
59. Utrero-Rico A, González-Cuadrado C, Chivite-Lacaba M, Cabrera-Marante O, Laguna-Goya R, Almendro-Vazquez P, et al. Alterations in circulating monocytes predict COVID-19 severity and include chromatin modifications still detectable six months after recovery. Biomedicines. (2021) 9. doi: 10.3390/biomedicines9091253
60. Rynikova M, Adamkova P, Hradicka P, Stofilova J, Harvanova D, Matejova J, et al. Transcriptomic analysis of macrophage polarization protocols: vitamin D3 or IL-4 and IL-13 do not polarize THP-1 monocytes into reliable M2 macrophages. Biomedicines. (2023) 11(2):608. doi: 10.3390/biomedicines11020608
61. Wigerblad G, Cao Q, Brooks S, Naz F, Gadkari M, Jiang K, et al. Single-cell analysis reveals the range of transcriptional states of circulating human neutrophils. J Immunol. (2022) 209:772–82. doi: 10.4049/jimmunol.2200154
62. Meng Y, Lin S, Niu K, Ma Z, Lin H, Fan H. Vimentin affects inflammation and neutrophil recruitment in airway epithelium during Streptococcus suis serotype 2 infection. Vet Res. (2023) 54:7. doi: 10.1186/s13567-023-01135-3
63. McColl SR, Paquin R, Ménard C, Beaulieu AD. Human neutrophils produce high levels of the interleukin 1 receptor antagonist in response to granulocyte/macrophage colony-stimulating factor and tumor necrosis factor alpha. J Exp Med. (1992) 176:593–8. doi: 10.1084/jem.176.2.593
64. Kunisada M, Yogianti F, Sakumi K, Ono R, Nakabeppu Y, Nishigori C. Increased expression of versican in the inflammatory response to UVB- and reactive oxygen species-induced skin tumorigenesis. Am J Pathol. (2011) 179:3056–65. doi: 10.1016/j.ajpath.2011.08.042
65. Xie X, Shi Q, Wu P, Zhang X, Kambara H, Su J, et al. Single-cell transcriptome profiling reveals neutrophil heterogeneity in homeostasis and infection. Nat Immunol. (2020) 21:1119–33. doi: 10.1038/s41590-020-0736-z
66. Glennon-Alty L, Moots RJ, Edwards SW, Wright HL. Type I interferon regulates cytokine-delayed neutrophil apoptosis, reactive oxygen species production and chemokine expression. Clin Exp Immunol. (2021) 203:151–9. doi: 10.1111/cei.13525
67. Leylek R, Idoyaga J. Chapter Four - The versatile plasmacytoid dendritic cell: Function, heterogeneity, and plasticity. In: Lhuillier C, Galluzzi L, editors. International Review of Cell and Molecular Biology. New York: Academic Press (2019). p. 177–211. Available at: https://www.sciencedirect.com/science/article/pii/S1937644819300991.
68. Varricchi G, de Paulis A, Marone G, Galli SJ. Future needs in mast cell biology. Int J Mol Sci. (2019) 20(18):4397. doi: 10.3390/ijms20184397
69. Nascimento C, Ferreira F. Tumor microenvironment of human breast cancer, and feline mammary carcinoma as a potential study model. Biochim Biophys Acta Rev Cancer. (2021) 1876:188587. doi: 10.1016/j.bbcan.2021.188587
70. Andrade LM, Sousa ECM, Melo ACS, Pardavil GSS, Silva RA, Valente KF, et al. Gene expression of cytokines in the blood of domestics cats infected and not infected by Feline alphaherpesvirus. doi: 10.4238/gmr19190
71. Smithberg SR, Fogle JE, Mexas AM, Reckling SK, Lankford SM, Tompkins MB, et al. In vivo depletion of CD4+CD25+ regulatory T cells in cats. J Immunol Methods. (2008) 329:81–91. doi: 10.1016/j.jim.2007.09.015
72. Rigamonti A, Castagna A, Viatore M, Colombo FS, Terzoli S, Peano C, et al. Distinct responses of newly identified monocyte subsets to advanced gastrointestinal cancer and COVID-19. Front Immunol. (2022) 13:967737. doi: 10.3389/fimmu.2022.967737
73. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. (2015) 33:495–502. doi: 10.1038/nbt.3192
74. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. (2019) 20:296. doi: 10.1186/s13059-019-1874-1
75. Zappia L, Oshlack A. Clustering trees: a visualization for evaluating clusterings at multiple resolutions. Gigascience. (2018) 7. doi: 10.1093/gigascience/giy083
76. Ahlmann-Eltze C, Huber W. glmGamPoi: fitting Gamma-Poisson generalized linear models on single cell count data. Bioinformatics. (2021) 36:5701–2. doi: 10.1093/bioinformatics/btaa1009
77. Yang S, Corbett SE, Koga Y, Wang Z, Johnson WE, Yajima M, et al. Decontamination of ambient RNA in single-cell RNA-seq with DecontX. Genome Biol. (2020) 21:57. doi: 10.1186/s13059-020-1950-6
78. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. (2019) 8:329–37. doi: 10.1016/j.cels.2019.03.003
79. Song Y, Miao Z, Brazma A, Papatheodorou I. Benchmarking strategies for cross-species integration of single-cell RNA sequencing data. Nat Commun. (2023) 14:6495. doi: 10.1038/s41467-023-41855-w
80. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141
81. Wickham H. ggplot2. New York: Springer International Publishing (2023). Available at: https://link.springer.com/book/10.1007/978-3-319-24277-4.
82. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. (2014) 32:381–6. doi: 10.1038/nbt.2859
83. R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical ComputingVienna, Austria (2023). Available at: https://www.R-project.org/.
84. Gu Z, Gu L, Eils R, Schlesner M, Brors B. circlize Implements and enhances circular visualization in R. Bioinformatics. (2014) 30:2811–2. doi: 10.1093/bioinformatics/btu393
85. Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, et al. Welcome to the tidyverse. J Open Source Software. (2019) 4:1686. doi: 10.21105/joss.01686
86. Chang H, Ashlock DA, Graether SP, Keller SM. Anchor Clustering for million-scale immune repertoire sequencing data. BMC Bioinf. (2024) 25:42. doi: 10.1186/s12859-024-05659-z
87. Radtanakatikanon A, Moore PF, Keller SM, Vernau W. Novel clonality assays for T cell lymphoma in cats targeting the T cell receptor beta, T cell receptor delta, and T cell receptor gamma loci. J Vet Intern Med. (2021) 35:2865–75. doi: 10.1111/jvim.16288
Keywords: feline, T cells, single cell RNA-sequencing (scRNA-seq), T-cell receptor repertoire, cross species analysis, myeloid Cells, V(D)J
Citation: Ramarapu R, Wulcan JM, Chang H, Moore PF, Vernau W and Keller SM (2024) Single cell RNA-sequencing of feline peripheral immune cells with V(D)J repertoire and cross species analysis of T lymphocytes. Front. Immunol. 15:1438004. doi: 10.3389/fimmu.2024.1438004
Received: 24 May 2024; Accepted: 23 September 2024;
Published: 15 November 2024.
Edited by:
Robert David Miller, University of New Mexico, United StatesReviewed by:
Kimberly Morrissey, University of New Mexico Health Sciences Center, United StatesDaron Standley, Osaka University, Japan
Copyright © 2024 Ramarapu, Wulcan, Chang, Moore, Vernau and Keller. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Stefan M. Keller, smkeller@ucdavis.edu