Comparison of single-cell landscapes among three amniotes reveals conservation and innovation in avian immune systems

doi:10.21203/rs.3.rs-5296203/v1

Download PDF

Article

Comparison of single-cell landscapes among three amniotes reveals conservation and innovation in avian immune systems

https://doi.org/10.21203/rs.3.rs-5296203/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Molecular characterization of chicken cells is essential for understanding avian physiology and vertebrate evolution, yet an organism-wide single-cell atlas in chicken is still lacking. Here we describe a comprehensive reference atlas of the chicken, encompassing 1.57 million cells across 149 cell types, along with a spatial transcriptomic map of the embryo. By integrating it with 1.2 million single cells from humans and turtles, we systematically explored the evolutionary rates of various cell types, particularly immune cells. The rapid evolution of chicken cells was generally characterized by changes in their gene regulatory networks and subsequent functional adaptations. In chicken, follicular dendritic cells emerge at the early development stage and exhibit myeloid rather than stromal origins, unlike in mammals. These cells share a regulatory network with mononuclear phagocytes and promote B cell proliferation and migration in the chicken-specific bursa of Fabricius. The observed variation in subtypes and proportions of γδ T cells across the three species reflected the evolution of pathogen recognition and signaling mechanisms among amniotes. Overall, our study provides an invaluable resource to study chicken cell biology and avian evolution and shines light on the evolutionary and cellular characteristics of immune cells across amniotes.

Biological sciences/Genetics/Genomics/Transcriptomics

Biological sciences/Immunology/Immunogenetics

Biological sciences/Genetics/Gene expression

Birds and mammals diverged approximately 310 million years ago¹, and birds have since developed several cells and tissues with specialized functions and physiological characteristics, including the gizzard, nucleated erythrocytes and the bursa of Fabricius. As the most species-rich class of amniotes and the closest phylogenetic relatives to mammals, birds represent a critical branch in the study of vertebrate evolution. Several bird species, particularly chicken (Gallus gallus), have been widely used as models for studying avian biology and evolution. In addition, chicken is well-recognized as a model for human biology and diseases, such as neuroscience, immunology, and development, due to the ease of accessing and manipulating living embryos². For example, HOXD11 has been identified as essential for limb bud development using chick embryos³. Moreover, birds are the major host of zoonotic pathogens such as avian influenza virus⁴ and Salmonella⁵, which are transmitted globally and deleterious to public health. A comprehensive characterization of avian cells, particularly immune cells, and their unique functions is thus a prerequisite for developing effective and sustainable strategies to control these infectious diseases.

Studies in birds have revealed notable differences in cellular function and composition compared to other animals, such as the ability of avian erythrocytes to recognize pathogens⁶, the presence of bursa-specific secretory dendritic cells (BSDCs)⁷, and a higher frequency of γδ T cells in peripheral blood⁸. Yet, investigating functions in avian cells using conventional cell-sorting approaches remains challenging due to their distinct cell receptors and the lack of reliable antibodies against them. Single-cell RNA sequencing(scRNA-seq) technology provided an unprecedented opportunity for characterizing cell composition, interactions, heterogeneity, and functions at the whole-transcriptome resolution without the need for species-specific antibodies. Single-cell atlases have been established in many species, including humans⁹, Macaca fascicularis¹⁰, and pigs¹¹. However, previous single-cell studies in chickens were severely limited by cell numbers, tissue types, and developmental stages^12–15. For example, scRNA-seq has been used to characterize the leukocyte profile in chicken peripheral blood¹², but circulating immune cells represent only a subset of the entire immune cell landscape and have not been compared with other amniotes. A comprehensive chicken cell atlas is thus essential to assess the avian cellular evolution and conservation by comparing it with those of other amniotes.

Here, we present a spatial transcriptomic map of chicken embryos and a comprehensive single-cell transcriptomic atlas using scRNA-seq, comprising over 1,576,581 cells from 36 tissues of chicken embryos and adults. To further explore the cellular evolution of amniotes and the unique characteristics of birds, we performed scRNA-seq on 232,761 cells from 14 tissues of adult turtles (Trachemys scripta elegans). By integrating publicly available 968,068 human cells from 32 tissues, we systematically studied the transcriptional evolution and conservation at single cell resolution across the three amniotes (Fig. 1A). Additionally, to validate the observed evolution of immune cells is avian-specific rather than chicken-specific, we newly generated single-cell data of 21,798 cells from three immune tissues of ducks (Anas platyrhynchos, i.e., spleen, bursa of Fabricius and intestine) to investigate the characteristics of the avian immune system. Altogether, these single-cell datasets generated from chickens, ducks and turtles provide a foundational resource for studying avian biology and vertebrate evolution at single cell resolution.

Construction of single-cell and spatial transcriptome atlas for three amniotes

To construct the chicken cell landscape, we performed sc/snRNA-seq across 10 and 35 tissues from three embryo and three adult animals, respectively, including immune-related tissues such as the thymus, caecal tonsil, bursa of Fabricius, spleen, and peripheral blood mononuclear cells (PBMCs) (Fig. 1B). After filtering out low-quality cells, we obtained a total of 1,576,581 single cells. The average number of cells and nuclei obtained from each tissue ranged from 5,084 cells (yolk sac) to 100,524 cells (gizzard), representing 0.32% and 6.4% of the total cells, respectively (Extended Data Fig. 1A). After dimensionality reduction and clustering of all the cells, 149 cell types were annotated based on the expression of cell canonical marker genes (Extended Data Fig. 1B and Supplementary Table 1). On average, 19 distinct cell types were identified per tissue, ranging from 6 in the testis to 39 in the spleen (Extended Data Fig. 1A). By comparing our chicken heart datasets with those published previously¹⁶, we found similar cell types between these datasets despite differences in developmental stages (Extended Data Fig. 1C-E). To study the spatial localization of cell populations, we generated spatial transcriptomics from the 9-day-old chicken embryo when the major organs appeared¹⁷. In total, we retrieved transcriptomic information for 76,154 bins, with an average of 19,281 captured genes per bin (Extended Data Fig. 2A). Unsupervised spatially constrained clustering of these bins showed transcriptomic configurations matching the localization of primary tissues (e.g., heart, lung and liver) (Fig. 1C, Extended Data Fig. 2B). The primary cell types exhibited a good match with spatial distribution in tissue anatomy (Fig. 1D, Extended Data Fig. 2C). For example, in the gizzard, pit cells were localized to the center, whereas the two types of smooth muscle cells and fibroblasts consistently mapped to the periphery (Extended Data Fig. 3A-C).

For the turtle cell landscape, we generated a total of 232,761 cells from 14 adult tissues, representing 59 distinct cell types (Fig. 1B, Extended Data Fig. 4A, B, and Supplementary Table 2). On average, 20 distinct cell types were identified per tissue, ranging from 9 in the testis to 42 in the lung (Extended Data Fig. 4A). Ionocytes, which specifically express FOXI1 and are important for ion transport and fluid pH regulation¹⁸, were found in the turtle lung and constituted 1.02% of all lung cells (Extended Data Fig. 4C, D). Out of 33,289 cells in the chicken lung, no ionocytes were identified. Whereas ionocytes were identified recently and comprised a low proportion of epithelial cells (0.01%) in human lung(~ 75,000 cells)¹⁹, indicating changes in cell type composition of amniotes from aquatic to terrestrial environments.

For human single cell data, we downloaded human 968,068 sc/snRNA-seq cells from 32 adult tissues that were nearly equivalent to those in chickens (Fig. 1B). We then analyzed and integrated them with our turtle and chicken data using 12,546 one-to-one orthologous genes. To address potential disparities in cell numbers among different cell types and discrepancies in labeling across subtypes between species, we excluded proliferative cells and consolidated subtypes within the same cell type. We performed downsampling for each cell type during the integration. The same cell types from the three species clustered together (Fig. 1E, Extended Data Fig. 5), suggesting the global expression patterns were highly similar among turtles, chickens and humans.

Comparison of cell landscapes among chickens, turtles and humans

Approximately 75% of orthologous genes were expressed in the same cell types across chickens, turtles, and humans (A gene was defined as expressed if it was detected in more than 25% of cells within a specific cell type and its expression level exceeded the third quartile for that cell type), while some genes were expressed exclusively in one species (Fig. 2A and Extended Data Fig. 6). In general, expression levels of orthologous genes of matched cell types in chickens were more similar to those in turtles compared to humans (Fig. 2B). Some cell types, such as neurons and erythrocytes, exhibited rapid evolution among amniotes. For erythrocytes, genes related to interferon signaling and antigen presentation exhibited species-specific expression in chickens (Fig. 2C, D), but were either lowly expressed or not expressed in turtles and humans, respectively. This was consistent with previous reports that nucleated erythrocytes in birds participate in immune responses and react to microbes and pathogens²⁰. Erythrocytes are also nucleated in turtles. However, the changes in erythrocyte function from turtles to chickens may reflect the enhancement of cellular immune functions during the transition from ectothermy to endothermy.

For cell types shared between chickens and humans, several adrenal cell types (e.g., chromaffin cell, capsular cells, and zona glomerulosa cell) showed lower correlations compared to other shared cell types (Fig. 2E and Extended Data Fig. 6A). Chromaffin cells are primarily associated with the synthesis and secretion of catecholamines²². In chickens, 1,468 genes exhibited up-regulated expression in chromaffin cells compared to humans (Supplementary Table 3). Many of these genes participated in the synthesis and secretion of catecholamines (SYT1 and SNAP25)^22,23 and the response to hormone (CHRM3 and GHR)^24,25 (Fig. 2E). Functional enrichment analysis revealed these genes were significantly enriched in receptor tyrosine kinases, hormone secretion and response (Fig. 2E). For zona glomerulosa cells, 1,471 genes were upregulated in chickens compared to humans (Supplementary Table 3). These genes, such as ACACA, FASN, HSD11B2, and HSD3B1, were enriched in response to hormones, lipid biosynthetic process, and metabolism of steroids^26–29(Fig. 2E). This indicates that chicken chromaffin cells and zona glomerulosa cells may possess elevated capacities for catecholamine and aldosterone production and secretion. Adrenal hormones play a critical role in regulating metabolism and electrolyte balance, including the elevation of blood glucose levels³⁰. The increase in hormone secretion levels may represent a gene regulation required during powered flight, an energetically demanding transport form. In addition to these two cell types, capsular cells up-regulated 1,392 genes in chickens compared to humans (Supplementary Table 3). These up-regulated genes were enriched in gland development, steroid hormone-mediated pathways, and NOTCH and WNT pathways (Supplementary Table 3). Capsular cells for protective outer layer of the adrenal gland, with NOTCH and WNT pathways supporting their proliferation and self-renewal^31,32. In addition to significant differences in cellular transcription, substantial differences in adrenal cell type composition were also observed between chickens and humans. In humans, the adrenal cortex is composed of three functional cell types: zona glomerulosa cells (ZG), zona fasciculata cells (ZF), and zona reticularis cells (ZR), which produce specific aldosterone, cortisol, and adrenal androgens, respectively³³. In contrast, only two subpopulations of adrenal cortex cells were identified in chickens (Fig. 2F). One cluster was highly enriched for the ZG marker genes HSD3B1 and AGTR1, and another cluster enriched for the ZF and ZR marker genes CYP11A1 AKR1B1 and CYB5A (Fig. 2G).

Cellular divergence driven by evolutionary lability in gene expression localization

Mononuclear phagocytes and lymphocytes are critical immune cell populations. Compared with humans, peculiarities in the chicken immune cells have been identified, such as the unique population of bursal secretory dendritic cells (BSDCs) in the bursa of Fabricius and the higher proportion of peripheral γδ T cells³³. In total, we obtained 56,286, 40,420 and 264,141 immune cells in the adult chicken, turtle and human, respectively. This large dataset allowed us to explore the evolution of immune cells in amniotes. In total, 30, 17, and 31 mononuclear phagocyte and lymphocyte cell types were identified in chickens, turtles and humans, respectively. PDCs were present in both chickens and humans, while DC2s, migratory DCs, and MAIT cells were found exclusively in humans. In addition, some shared immune cell types, such as B cells and pDCs, exhibited distinct tissue distribution patterns. Besides the discrepancy in immune cell composition and distribution, we next quantified the similarity among the average transcriptomes of the immune cell types. We first used pairwise unsupervised MetaNeighbor analysis and the mean AUROC score to quantify the similarity between cell-type pairs³⁴. Although human FDCs are stromal-derived cells rather than myeloid-derived cells, we included them in the following analysis to assess their similarities to chicken FDCs in gene expression. Cell-type dendrogram clustered immune cells into three major categories: mononuclear phagocytes, B lymphocytes, and T/innate lymphocytes (Fig. 3A). While most cell types were arranged in accordance with the categories, some immune cell subtypes exhibited deviations, especially plasmacytoid dendritic cells(pDCs) (Fig. 3A). Human pDCs clustered with B cells, whereas chicken pDCs clustered with mononuclear phagocytes. Further investigation revealed that 1,695 genes exhibited significant expression divergence in pDCs between the two species, with 898 showing a higher expression in chickens and 797 in humans (Supplementary Table 4). These genes with explicit expression in chicken pDCs were significantly enriched in viral infection, protein tyrosine kinase activity, platelet activation and leukocyte migration (Supplementary Table 4). Evolutionary changes in gene expression suggested that many of these genes showed marked changes in cell-type localization of expression between chicken and humans (Fig. 3B and C, respectively). For example, genes related to lipid and lipoprotein transport, and those involved in scavenging free heme for iron recycling, were expressed in chicken pDCs, but in human erythrophagocytic macrophages (Fig. 3B, C). Additionally, genes related to leukocyte migration and adhesion were expressed in chicken pDCs, but in human migratory dendritic cells(migDCs) (Fig. 3B, C). These findings suggested that pDCs in chickens and humans underwent relatively rapid cellular divergence, driven partially by change in gene expression modules.

To systematically identify evolutionary changes in the expression of orthologous genes in immune cells of chickens and humans, we referred to a previous method and divided expression patterns into four types³⁵: Type 0 ('conserved'), Type 1('expression gain/loss'), Type 2 ('expression expansion/contraction') and Type 3('expression switch'). Through comparing gene expression of the same cell types between chickens and humans, Type 0 denotes genes that are expressed in both species, whereas type 1 denotes genes that are expressed in only one. Type 2 changes involve the gain or loss of expression in additional cell types during evolution. Type 3 ('expression switch') changes involve a switch in expression from one cell type to another. Only 5.05% of expressed genes showed fully conserved patterns between species. Most genes underwent expression gain/loss (type 1, 46.8%) or expansion/contraction (type 2,40.8%) (Supplementary Table 5). Given the divergence in pDC transcription observed above, we further focused on evolutionary changes in gene expression within mononuclear phagocytes. Only 9.49% of expressed genes showed fully conserved patterns (Fig. 3D and Supplementary Table 6), suggesting expression patterns of nearly all genes are evolutionarily labile. For example, Type 0 gene CR2 was only expressed in FDCs in chickens and humans (Fig. 3E). In contrast, Type 1 gene CXCL14 was expressed only in FDCs in humans (Fig. 3E). For type 2 and type 3 genes, two important examples were ITGB2 and LYZ. ITGB2, a known regulator of phagocytosis³⁶, was mainly expressed in chicken FDCs and DC1s but in monocytes, macrophages, and DC1s in humans. LYZ encoding the antibacterial enzyme was selectively expressed in pDCs and FDCs in chickens, but in monocytes, macrophages, and DC1s in humans (Fig. 3E).

Avian FDC fulfill cellular migration and development for both B cell progenitors and germinal centre B cells

In the chicken bursa of Fabricius, bursal secretory dendritic cells (BSDCs) are thought to contribute to the microenvironment necessary for B-lymphocyte differentiation³⁷. However, their specific characteristics and role within follicle buds remain poorly understood. Here we analyzed 34,861 single cells from the embryonic bursa of Fabricius. One distinct cell cluster demonstrated a high and specific expression of BSDCs marker genes (CSF1R and CD74)³⁸, but not macrophage-specific markers C1QA and C1QC (Fig. 4A, B). Spatial transcriptomics on the embryonic bursa of Fabricius further revealed that this cluster was scattered within the follicle and colocalized with B cells (Fig. 4C, D and Extended Data Fig. 7A), indicating that this cluster is BSDCs. We then compared the transcriptional patterns of BSDC with those of other chicken mononuclear phagocytes. BSDCs in chickens specifically and highly expressed ZNF366, LAMP3, SPIC, MAFB and LYZ (Fig. 4E), which are involved in the maturation and activation of mammalian DCs^39,40, the differentiation of macrophage and antibacterial activity^41,42. Among them, SPIC exhibited significant regulatory effects (Extended Data Fig. 7B). Up-regulated genes in BSDCs were significantly enriched in phagocytosis, adhesion, and positive regulation of leukocyte cell-cell adhesion (Extended Data Fig. 7C). Interestingly, they also expressed mammalian FDC marker genes CXCL13, CR2 and TNFSF13B and interacted with B cell progenitors, partly mediated through CXCL13 and TNFSF13B (Fig. 4F). CXCL13 was the main chemokine secreted in BSDCs, and its receptor CXCR5 was expressed by B cell progenitors (Fig. 4F). In mammals, CXCL13/CXCR5 interaction from FDCs to mature B cells promotes the migration of the latter^43,44. Parallel to this process, CXCL13/CXCR5 interaction from BSDCs to B cell progenitors facilitates the migration of these progenitors into the follicles of the bursa. In the bursa, only B cells not other immune cells, can migrate into the follicles⁴⁵, likely attributed to the unique interactions between B cells and BSDCs. Besides, BAFF released by BSDCs interacted with BAFF-R in B cell progenitors (Fig. 4F). This signaling promotes B cell proliferation and aids in maintaining the B cell population by preventing apoptosis⁴⁶. The above results indicated that BSDCs play an indispensable role in B cell development, primarily through the recruitment of B cells and regulatory function. Notably, in chickens, B cell lymphopoiesis is confined to the hematopoietic organs (such as the spleen) and the bursa of Fabricius, which was in line with the distribution of BSDC cells (Extended Data Fig. 8). By contrast, human B cell progenitors mainly develop in the bone marrow⁴⁷. The indispensable function and high tissue-specific distribution of BSDC partially explain why the bursa is the sole site for B cell development in chickens.

Given that BSDCs exhibit high expression of FDC marker genes, we further explored the relationship between these two cell populations. In adult chickens, FDCs were identified using well-established marker genes CXCL13, CR2, and TNFSF13B. The distribution of FDCs aligned with that of GC B cells, predominantly found in the caecal tonsil, cecum, and bursa (Extended Data Fig. 8). The spatial transcriptomics data of the caecal tonsil showed that FDCs and GC B cells shared spatial neighbourhoods, indicating a preferential interaction between FDCs and B cells (Fig. 4H). This interaction, mediated by CXCL13-CXCR5, VCAM1_ITGA4(VLA-4 subunit alpha) and TNFSF13B-TNFRSF13B signaling, is the same as the mechanism between human FDC and GC B cells (Extended Data Fig. 7H). These findings also support the accuracy of the annotation for chicken FDCs. When examining all mononuclear phagocytes, FDCs and BSDCs clustered together (Extended Data Fig. 7D-F). Correlation analysis revealed that the transcriptional patterns of the two cell types were highly similar (R = 0.92) (Fig. 4G), indicating that BSDCs and FDCs are the same cell types at different developmental stages. Up-regulated genes in embryos were mainly enriched in protein transferring and processing, as well as histones/ribosome binding (Extended Data Fig. 7G). Conversely, up-regulated genes in adults showed enrichment for immune response and antigen presentation (Extended Data Fig. 7G). This transcriptional transition accompanies the maturation of immune functions analogous to the developmental process of mononuclear phagocytes in mammals⁴⁷. Overall, chicken FDCs fulfill a ‘double-duty’ role in cellular migration and development for both embryonic B cell progenitors and adult germinal center B cells through similar mechanisms (Fig. 4I).

Human FDCs are differentiated from perivascular precursors of stromal cells, whereas chicken FDCs are differentiated from hematopoietic stem cells⁷. Pseudotime analysis also confirmed the myeloid origin of chicken FDCs (Extended Data Fig. 9A and Supplementary Table 7). To gain insight into the relationship between FDCs and other mononuclear phagocytes of amniotes, we integrated a total of 32,379 mononuclear phagocytes and FDCs from turtles (8,488; 0), chickens (11,517; 959) and humans (12,374; 203). The integrated cell map showed that monocytes, macrophages, DC1s and pDCs were distributed in the same cell clusters across the species (Fig. 4J). By contrast, FDCs in chickens and humans were distributed in entirely different cell clusters. Human FDCs were concentrated in cluster 3, whereas chicken FDCs were in cluster 10 (Fig. 4J). Moreover, chicken FDCs were closer to other mononuclear phagocytes, rather than human FDCs (Fig. 4J). Pairwise cluster correlation analysis revealed that FDCs in chicken were more similar to myeloid cells than to human FDCs (Extended Data Fig. 9B). The above results indicated that, although FDCs in chicken and human expressed the same marker genes (CXCL13 and CR2), their ontogeny and core regulatory complexes (CoRCs) were different. Human FDCs expressed marker genes TMEM119 and SOX9⁴⁸,⁴⁹(Fig. 4J), whereas chicken FDC expressed ZNF366, SPIC and MAFA (Fig. 4K). In addition, by examining 17,383 cells from the duck spleen and bursa of Fabricius, we found that duck FDCs exhibited greater similarity to human monocytes than to FDCs (Extended Data Fig. 9C, D). In summary, the presence of distinct CoRCs enabled us to identify FDCs in both birds and mammals, which appear to have evolved similar functions through evolutionary convergence.

The evolution of GC B cells and affinity maturation of antibodies in amniotes

Germinal centers (GCs), present in birds and mammals, are key histological structures for generating high-affinity B cells. Reptiles lack visible GCs in their spleen and show limited affinity maturation and secondary antibody response⁵⁰. However, sharks exhibit GC-like B cells and affinity maturation in splenic follicles⁵¹, prompting questions about the evolution of GCs and affinity maturation in vertebrates. Here we conducted a deep analysis and comparison on B cells in turtles (n = 13,686)) and chickens (n = 38,810). In adult chickens, we detected naive B, memory B, plasma B, and two populations of germinal center B cells (GC_B(I), GC_B(II)), consistent with mature and antigen-experienced B cell types in humans (Fig. 5A). Naive B cells and memory B cells were characterized by differential expression of KLF2, GPR183, BCL2 and HHEX^52–55(Extended Data Fig. 10A). In turtle, only four cell types or states were identified (Fig. 5A), including naive B, activated B, cycling B and plasma B cells. Naive B cells and activated B cells were distinguished by differential expression of KLF2, IL16 and BCL6(Extended Data Fig. 10A). Cell type annotation indicated the absence of GC_B(II) cells in turtles. To further validate this at the transcriptional level, we compared B cells of turtles and chickens in detail (turtle:13,686; chicken:12,000, downsampled from the population of 38,810 B cells to make it comparable to the turtle dataset), due to the evolutionary proximity between them⁵⁶. Pairwise correlation analysis confirmed that among all B cells, chicken GC_B(II) cells exhibited the low similarity with turtle B cells and lacked a good counterpart in turtles. Integrative analysis also indicated that B cells from chicken co-clustered with their counterparts in turtle, and GC(I)_B of chicken co-clustered with cycling B cells of turtle (Extended Data Fig. 10B, C). Only GC(II)_B, concentrated in cluster 1, did not have a good counterpart in turtles. We also examined the expression of genes related to germinal center response in amniotes, including AICDA for somatic hypermutation and class-switch recombination⁵⁷. The results showed that turtles exhibited lower levels of AICDA and BCL6 expression compared to those in chickens and humans (Fig. 5C). This supports the Ag-receptor diversity of B cells in reptiles, but they may lack robust affinity maturation^58,59. Interestingly, in amphibian Xenopus⁶⁰, the expression of AICDA was comparable with that in turtles. However, MKI67, a reliable marker of proliferating cells in Xenopus, was scarcely expressed, and no cluster of MKI67-enriched proliferating B cells was observed (Fig. 5C). Analogs of GC B cells in sharks suggested the evolutionary foundation of GCs dates back to the jawed vertebrate ancestor. During animal evolution, the number of B cell types has changed. GC_B(II) cell types are lost in some amphibian and reptilian clades, such as Xenopus and turtle.

Although the primary B cell types in chickens and humans were similar in adults, the proportion of subsets showed significant differences between these two species. Compared with humans, the proportion of GC B cells was remarkably increased in chickens, accompanied by a decrease in memory cells, probably reflecting high germinal center activity (Fig. 5D). The significant difference in proportions prompted us to explore the differences in the development of B cells between chickens and humans. We compared gene expression levels in chicken cell types with those of their orthologous genes in the corresponding human cell types(chicken 12,000 B cells and human 9,022 B cells). GC B(II) cells in the two species consistently exhibited low correlation, irrespective of the number of variable genes used (Fig. 5E). Integrated analysis of chicken and human B cells also validated the lower similarity of GC(II)_B (Extended Data Fig. 10D, E). Gene expression analysis of GC(II)_B between chicken and human revealed 345 genes with species-biased expression patterns, including 121 up-regulated genes and 222 down-regulated genes in chicken (Fig. 5F and Supplementary Table 8). Species-biased genes were enriched in lymphocyte activation and cell cycle categories, such as cell responses to stress, B cell receptor signaling pathway, and positive regulation of programmed cell death (Fig. 5G). For example, MITF is a negative regulator of BCR signaling⁶¹. The voltage-gated proton channel HVCN1 promotes the production of ROS, which augments the proliferation of activated B cells and delays plasma cell differentiation^62–65. In the germinal centre, B cells interact with FDCs and T follicular helper cells(Tfhs), which favours the survival of higher affinity B cells and forces others to undergo apoptosis by neglect⁴⁸. Expression differences in genes related to BCR signal transduction and affinity selection pressures suggested potential differences in the production of high-affinity antibodies between chickens and humans. Cell type-specific CoRCs are the driver of cell type identity⁶⁶. Transcriptional regulatory network analysis for GC(II)_B in chickens and humans showed that they shared transcription factors and regulatory mechanisms for GC B cell activation and survival, such as BCL6 and PAX5 (Fig. 5H). Additionally, vital regulators of BCR signaling were unique in chickens, and genes encoding these regulators were highly expressed in chicken B cells. To validate whether these unique transcription patterns also exist in other birds, we compared the gene expression in ducks and humans' GC B(II) cells. Results showed that genes with high expression in chickens were also highly expressed in ducks (Fig. 5I). This suggests that GC B(II) cells have undergone rapid molecular evolution in birds compared to other mature B cell subsets.

Divergent evolution of the γδ T lineage in amniotes

Previous studies have revealed the T cells and innate lymphocyte subpopulations in chicken PBMCs using scRNA-seq⁶⁷. Canonical marker genes of these cells are conserved between chickens and humans. In our comprehensive immune cells profiling, we identified 33 transcriptionally distinct cell populations in both embryonic and adult chickens (Fig. 6A). Most cell types were consistent with those found in mammals, except for two unique clusters. One cell cluster, named CTSG + immune cells, was highly enriched for the lymphocyte marker genes such as CD3D and CD4, as well as genes preferentially expressed in mammalian granulocytes or mast cells(CTSG, NDST2, HDC, and CSF2RB) (Fig. 6B)⁶⁸. This specific gene combinatorial code diverged from that of known conventional human immune cells. Up-regulated genes in CTSG + immune cells were enriched in leukocyte activation and differentiation (Fig. 6C). Transcriptional pairwise cluster correlations between CTSG + immune cells and other immune cells revealed that it showed the strongest mutual correlations with ILC2 and heterophil, the avian equivalents of neutrophils (Fig. 6D). Another cluster was enriched with genes involved in T cell activation, positive regulation of cytokine production and cytoskeleton regulation, but almost did not express CD4 and CD8A. This suggests that this cell population may be a non-canonical T cell subset with functions mediated through cytokines (Fig. 6E). Notably, this cell cluster was only localized to adult caecal tonsil and spleen, indicating high tissue specificity (Extended Data Fig. 8). Compared to chickens and humans, only 13 T cell types/states were identified in turtles (Fig. 6F), including major cell types such as ILC, NK, γδT and CD4+/CD8 + T cells.

To examine T cells and innate lymphocytes with conserved or innovative gene expression profiles in adult turtles, chickens, and humans, we integrated 16,998, 29,800, and 28,258 cells in turtles, chickens and humans, respectively. Clustering of the integrated data yielded 37 cell clusters (Fig. 6G-H). Generally, CD4+/CD8 + T cells and NK cells showed high similarity among amniotes (Extended Data Fig. 11A). For example, naïve T cells from these three species co-clustered and segregated into 4 clusters (C12, C37, C1, C35) (Fig. 6I-J, Extended Data Fig. 11A). They shared expression of several transcription factors, including KLF2, TCF7, SATB1 and FOXP1(Fig. 6J), confirming that the CoRCs of naïve T cells were conserved across amniotes. Similarly, we also observed Tregs from these three species co-clustering in the same neighborhoods (cluster 7,11) (Fig. 6I), which expressed CTLA4 at high levels (Fig. 6J). Conversely, several integrated clusters included cells from chickens only, indicating that these cell types have unique gene expression profiles. For example, cluster 32 was mainly enriched in chicken γδ T_1 cells(Fig. 6H). Consistent with this, the gene combinations specifically expressed in this cell type were not detected in either turtle or human (Fig. 6J). Moreover, other γδ T cells from these three species did not co-cluster either, and the effector genes were also different (Fig. 6H and Fig. 6K, Extended Data Fig. 11A). For example, γδ T cells in humans with a distinctive expression of the cytotoxic effector molecul GZMA, co-clustered with effector CD8 + T cells. In turtles, γδ T cells mainly expressed the cytotoxic effector molecule GZMH. In chicken, γδ T_2 cells were committed to producing bacteriostatic or lytic molecules, such as IL17A and GNLY, and they co-clustered with Th17 cells. Weighted gene correlation network analysis also revealed that chicken γδ T_2 cells and Th17 cells shared a module and associated genes (Extended Data Fig. 11B, C). To explore if the main subpopulation of γδ T cells in other birds was consistent with those in chickens, we examined immune cell types and transcription patterns in the duck intestine. The results showed that duck γδ T cells were also enriched in IL17(Extended Data Fig. 11D, E). Taken together, these results indicated functional specialization of γδ T cell subpopulations across different species, driven by the divergency of core regulatory programs.

Avian γδ T cell function and ontogeny

Previous studies have shown that the frequency of chicken γδ T cells in peripheral blood is relatively high (20–50%)⁸, compared with that in humans (less than 5%)⁶⁹. Therefore, chickens are referred to as “γδ T cell high” species, while humans belong to “γδ T cell low” species. Based on our data, proportions of γδ T cells in chickens increased by approximately 6.7% in all tissues compared to that in humans (Fig. 7A). When focusing on PBMC and intestine, γδ T cells in chickens increased by approximately 20% and 8% compared to that in humans, respectively (Fig. 7A). Apart from this, the integrated analysis above also indicated that γδ T cell subtypes and functions were significantly different between chicken and human. This substantive difference prompted us to explore the function and ontogeny of avian γδ T cells.

In adult chickens, there were two γδ T cell subsets in peripheral tissue. Chicken γδ T_1 cell up-regulated genes that encode γδ T effector programming transcription factors SOX13 and MAF, as well as cytokine receptors (IL20RA, IL9R), but did not show cytokine production (Fig. 7G). SOX13 and MAF are essential for the establishment of γδ T cell identity and the commitment of IL-17-producing γδ T cells^70,71, indicating that γδ T_1 cells were not a terminally differentiated population. The γδ T_2 subpopulation up-regulated genes encoding bactericidal molecules (GNLY) and cytokines (IL17A, IL17F, IL22) (Fig. 7G), which are involved in protective immunity against extracellular bacteria and fungi. Pseudotime analysis predicted a trajectory from immature γδ T cells in the thymus through γδ T_1 to γδ T_2(Fig. 7B, C and Supplementary Table 9). Additionally, γδ T_1 cells were present in the thymus, PBMC and other tissues. In contrast, γδ T_2 cells were not identified in the thymus or PBMC, but only in other peripheral tissues (Extended Data Fig. 8). This suggested that γδ T cells in chicken commit to effector cytokine production in peripheral organs, rather than in the thymus. In peripheral tissues, apart from the direct function of pathogen clearance, the cytokines released by γδ T cells can interplay with other immune cells⁷², epithelial cells and fibroblasts to exert an immunoregulatory effect. To gain insight into the potential function of γδ T cells, we examined spatial transcriptomics in the chicken caecal tonsil, where immune cells were concentrated. Spatial distribution showed γδ T cells localized close to B cells. Local hotspots of ligand-receptor pairs occurred at the interface and mediated the interaction between them. These interactions included CD40LG-CD40, ICOS-ICOSL and TNFSF8-TNFRSF8 interaction(Fig. 7D), which support B cell activation, growth and differentiation⁷³. In mammals, IL17-expressing Th17 cells can function as B-cell helpers by not only triggering B cell proliferation but also promoting class-switch recombination⁷⁴. Interestingly, Th1/Th17 were approximately 23% lower in chicken intestine compared with those in humans (Fig. 7A). Therefore, γδ T cells expressing IL-17 function similarly to Th17 cells to exert immune regulatory effects in chickens. Additionally, γδ T cells were also physically close to the enterocytes. IL-17 released by γδ T_2 cells was predicted to interact with IL-17 receptors in enterocytes (Fig. 7D). This interaction could induce the production of proinflammatory cytokines and chemokines, thereby recruiting lymphocytes⁷⁵. In summary, the high frequency of γδ T_1 cells in adult chicken peripheral tissues contributes to the effective supply of effector γδ T_2 cells, which possess powerful antibacterial properties and bridge innate and adaptive immunity (Fig. 7E).

It has been reported that mammalian γδ T cells exhibit heterogeneity across developmental stages and tissues⁷⁶. Based on this, we explored the heterogeneity of chicken γδ T cells. In addition to adult γδ T cells, we identified three distinct γδ T cell subsets in embryos, all of which displayed transcriptional profiles closely resembling those of adult γδ T cells (Fig. 7F, G). For example, embryonic γδ T_2 cells highly expressed chemokines and IL-17 signaling molecules. They shared transcription modules with adult γδ T_2 subpopulation (Fig. 7G, Extended Data Fig. 12A). Gene expression comparison between them revealed that embryonic γδ T_2 cells up-regulated interferon-related genes (IRF8, IFIH1, IFIT5), while adult γδ T_2 cells up-regulated genes related to antigen processing and presentation (GPR183, CD82, TNFRSF9) (Fig. 7H, Extended Data Fig. 12B). This indicated that the adaptive immune function of adult γδ T cells was gradually refined. To comprehensively understand the heterogeneity across tissues, we focused on the adult γδ T cells, due to their widespread organ distribution. In adult chickens, γδ T_1 cells showed a high tissue heterogeneity (Extended Data Fig. 12C). Thymic γδ T_1 cells expressed gene rearrangement and lineage differentiation-related genes (RAG1, RAG2, SOX13 and TARP) (Fig. 7I). PBMC γδ T_1 cells overexpressed interferon and antigen presentation related genes (BF1, IFI6 and IRF1) (Fig. 7I). Splenic and intestinal γδ T_1 cells had higher expression of chemotaxis genes. In contrast, γδ T_1 in non-immune organs had higher expression of genes associated with T cell activation (Fig. 7I). These findings suggested that the maturation and activation of γδ T_1 cell accompany their egress from the thymus to seed peripheral tissues, where the local microenvironment shapes them into populations with distinct effector functions. In contrast, γδ T_2 could be further subdivided into three subsets (Extended Data Fig. 12D), all of which are distributed simultaneously in peripheral tissues. Cluster 1showed high expression of genes associated with lymphocyte-mediated antibacterial activity (GNLY), cluster 2 γδ T cells exhibited abundant pro-inflammatory cytokines and chemokines, whereas cluster 3 expressed stemness-associated markers (e.g., CCR7, TCF7) (Extended Data Fig. 12E).

Our study provides a comprehensive single-cell dataset from multiple chicken organs, spanning > 1,500,000 single-cell profiles from 36 tissues and identifying 149 cell types. Using spatial transcriptomics, we delineated tissue organization and cellular communication networks. Besides, we profiled 232,761 cells from 59 cell types in adult turtles, allowing us to systematically compare cell types across amniotes by integrating our findings with published human single-cell transcriptome data. Notably, we detailed the immune cell landscapes to characterize avian immune features and to investigate the evolution of immunity across amniotes.

Amniote vertebrates (reptiles, birds and mammals) originated from a common ancestor about 310 million years ago¹. Our findings revealed high similarities in cell types between turtles, chickens and humans, underscoring the conservation of gene expression patterns across amniotes. While similarities in cell types across these species are evident, their evolutionary rates vary due to distinct evolutionary pressures. Previous studies have revealed the evolution of neurons and testicular cells in terms of anatomical structures, cell types, and molecular patterns, as shown by single-cell sequencing of brains and testis^77–79. Our systematic comparison of amniotic cell landscapes confirmed rapid evolution in certain cell types, such as Sertoli cells and neurons, while also identifying other rapidly evolving cell types, including erythrocytes and adrenal cortex cells. Sequence divergence in non-coding genome regions likely drives the emergence of species-specific traits⁸⁰. Future research should explore the evolution of gene regulatory programs to better understand how genetic divergence contributes to species-specific phenotypes.

Accumulating studies demonstrate that molecular conservation is foundational to innate and adaptive immunity across vertebrates⁸¹. Consistently, in this study, we identified conserved immune profiles among amniotes in terms of immune cell populations and core regulatory complex, including macrophages, CD4+/CD8 + T cells, NK and B cells. Interestingly, innate immunity, considered primitive during evolution, displays low evolutionary conservation. Innate immune cell types have diversified during the evolution of amniotes. This diversification may reflect species-specific differences in pathogen recognition and signaling mechanisms. Furthermore, avian immune cell lineages, such as pDCs, FDCs, and γδ T cells, exhibit distinct characteristics compared to their mammalian counterparts.

FDCs, for example, are essential for lymphoid follicles organization and the germinal center reaction. In humans, FDCs differentiate from stromal cells with the assistance of mature B cells^82,83. However, avian FDCs were of hematopoietic origin, enabling their presence during embryonic hematopoiesis without mature B cells support. This early presence may offer two primary advantages for birds. First, it can assist B cell progenitors in migrating to the bursa of Fabricius (Fig. 4). Second, it may enhance the efficiency of germinal center formation and antibody production in birds. After hatching, chicks are immediately exposed to environmental antigens, pathogens, and gut microbiota colonization, increasing their immune system demands. The GC reaction is critical for production, as FDCs retains antigen-antibody complexes via complement and Fc receptors, supporting GC B cells survival⁸⁴. Therefore, developing the follicular dendritic cell (FDC) network influences the germinal center responses. Previous research showed that the first germinal centers during the chicken ontogeny appear as early as the fourth day after hatching⁸⁵. By contrast, there were no FDCs within the lymphoid microarchitecture in 7-days-old mice⁸⁶ and primary follicles first appeared on day 12 in the mesenteric nodes of mice⁸⁷. Therefore, the early emergence of avian FDCs may be crucial for initiating antibody responses post-hatch.

Apart from FDCs, innate γδ T cells also exhibit more rapid divergency than other immune cells in amniotes, varying in both proportions and global gene transcription. The high frequency of γδ T cells in both chicken embryos and adult chickens suggests an important role in pathogens' defence. In chicken embryos or young chicks, adaptive immune functions are relatively underdeveloped, positioning γδ T cells as critical components of the host's defence. In adults, γδ T cells are widely distributed across tissues in a preactivated or primed state, serving as frontline defenders. For example, increased numbers of γδ T cells have been reported in chickens' peripheral blood and spleen shortly after exposure to Marek’s disease virus and Salmonella^88,89. Interestingly, similar to chickens γδ T cells constitute up to 10–20% of the T-lymphocyte population in peripheral blood and are distributed throughout epithelial surfaces in adult ruminants⁹⁰. These γδ T cells also have an important role in many infectious diseases. Compared with humans, these farm animal species are γδ-T-cell-high species. Contrary to αβ T cells, γδ T cells recognize a broad range of antigens without restriction by major histocompatibility complex molecules and are primed for rapid effector function⁹¹. This characteristic may be particularly beneficial in defending farm animals against diverse environmental pathogens.

In conclusion, our single-cell atlases of chicken and turtle, combined with spatial transcriptomic data, provide valuable resources for future research into avian morphological features, vertebrate evolution and zoonotic disease mechanisms. These findings lay the groundwork for further studies into the evolutionary pathways that have shaped immune responses across species.

Ethics statement

The experimental protocol was approved by the Institutional Animal Care and Use Committee of HUNAU (2024 − 159). All experimental procedures were conducted following the national and institutional guidelines for using experimental animals for research.

Collection of chicken, duck and turtle tissues

A total of three adult White Leghorn chickens (Gallus gallus, approximately 14 weeks) and three adult Xiangjia Black Phoenix chickens (approximately 13 weeks), three chicken embryos (17–21 days), three adult Peking ducks (Anas platyrhynchos, approximately 24 weeks) and three adult red-eared slider turtles (Trachemys scripta elegans, approximately 5 years) were obtained from the farm. Chicken, duck and turtle tissues were isolated, rinsed with PBS and minced into small pieces by mechanical dissociation. Next, the fresh samples prepared for scRNA-seq were transferred to MACS® tissue storage solution (Miltenyi Biotec Technology & Trading, #130-100-008) and stored at 4℃ until they were dissociated within 48 hours. Samples for snRNA-seq were transferred to cryogenic vials and then frozen and stored in liquid nitrogen until nuclear extraction was performed. Peripheral blood mononuclear cells from heparinized venous blood were isolated using mononuclear cell isolation kit (Solarbio, P8910) according to the protocol. Then cells from the blood were resuspended in freezing medium composed of 90% FBS and 10% DMSO and frozen using a freezing container in a -80°C freezer for 24 hours before being transferred to liquid nitrogen for long-term storage.

Single-nucleus/cell suspension preparation

Single-nucleus isolation was performed as previously described^10,11. In brief, frozen tissues were transferred to homogenizer with 1 ml lysis buffer consisting of 250 mM sucrose (Ambion), 10 mg/mL bovine serum albumin (Ambion), 5 mM MgCl2 (Ambion), 0.12 U/µL RNasin Plus (Promega, #N2115), 0.12 U/µL RNasein (Promega, #N2115) and 1× cOmplete Protease Inhibitor Cocktail (Roche, #11697498001). After two additional rounds of homogenization, the mixture was filtered through a 40-µm cell strainer and centrifuged at 500g for 5 min at 4°C. The pellets were resuspended in 1 ml of buffer B containing 320 mM sucrose, 10 mg ml-1 BSA, 3 mM CaCl2, 2 mM magnesium acetate, 0.1 mM EDTA, 10 mM Tris-HCl, 1 mM DTT, 1× cOmplete Protease Inhibitor Cocktail and 0.12 U µl-1 RNaseIn. After centrifugation as described above, nuclei were resuspended with cell resuspension buffer at a concentration of 1,000 nuclei per µl for library preparation.

For cell suspension preparation, the fresh samples were incubated for 1 h in 10 ml DS-LT buffer (0.2 mg/ml CaCl2, 5 µM MgCl2, 0.2% BSA and 0.2 mg/ml Liberase in HBSS) at 37°C. After this, the tissue digestion was stopped by adding FBS. Then, the suspension was filtrated through a 100 µm cell strainer and centrifuge. Cells from the spleen were obtained from fresh tissue by mechanical dissociation. Cells from PBMCs were obtained as described above. Samples were filtered through a 40 µm cell strainer and centrifuged. Pellets were resuspended in cell resuspension buffer at 1,000 cells per µl for library preparation.

Single-cell library construction and sequencing.

The DNBelab C Series Single-Cell Library Prep Set (MGI, 1000021082) or 10X Chromium system were employed for library preparation according to the previously established protocols^10,92. Briefly, single-nucleus/cell suspensions were used for droplet generation, emulsion breakage, bead collection, reverse transcription and cDNA amplification to generate barcoded libraries. Indexed libraries were constructed according to the manufacturer’s protocol. After the measurement of cDNA concentrations, libraries were sequenced on DNBSEQ-T7 or Illumina PE150 platform.

scRNA-seq and snRNA-seq data processing

After the filtering of raw sequencing reads, reads were aligned to the chicken genome (GCF_016699485.2), duck genome(GCA_008746955.1) and turtle genome (GCF_013100865.1). Three files related to genes, barcodes and the raw UMI count were generated by the DNBelab C Series scRNA analysis software (https://github.com/MGI-tech-bioinformatics/DNBelab_C_Series_scRNA-analysis software) or Cell Ranger software. Then, we created a Seurat object using the Seurat R package⁹³. For each sample, cells were retained for downstream analysis if the number of detected genes exceeded 500 but was below the 95th percentile threshold, and if the percentage of mitochondrial genes was less than 15%. After filtering, the raw count matrix was normalized and scaled. The top 2,000 most variable genes of each sample were identified using the “FindVariableFeatures” function. Harmony performed batch correction across replicates⁹⁴. Dimension reduction, clustering and visualization were performed with default parameters following the guidance. Only genes with an adjusted P-value < 0.05, an absolute log fold change (LogFC) > 0.25, and detected in at least 25% of the cells within a cluster were identified as marker genes. Finally, each cluster was annotated based on the expression of established marker genes.

Spatial transcriptomics sequencing and data analysis

For chicken embryo, spatial transcriptomics was performed using Stereo-seq platform according to the standard protocol⁹⁵. Tissue section was adhered to the Stereo-seq chip surface and fixed in methanol, followed by nucleic acid staining and imaging. After permeabilization, mRNAs captured by DNA nanoballs (DNBs) on the chip were reverse transcribed and amplified. DNBs were then loaded onto the patterned Nano arrays and sequenced using the MGI DNBSEQ-Tx sequencer. After raw data was processed, we treated the bin 100 as the analysis unit and performed unsupervised clustering. Data normalization, scaling, and bins clustering were processed using the R package Seurat⁹³. Tissue identities of clusters were annotated using tissue-specific expression genes. For the spleen and gizzard, we extract the corresponding areas and performed spatial reclustering. Specifically, we treated bin30 of the spleen and bin50 of the gizzard as the analysis units.

For bursa of Fabricius and caecal tonsil, spatial transcriptomics was carried out using the BMKMANU S1000 platform according to the protocol⁹⁶. Tissue sections of the chicken embryo bursa of Fabricius and adult chicken caecal tonsil were placed on sequence slides. The sections were then fixed with cold methanol, followed by H&E staining and imaging before permeabilization, according to the user guide for the BMKMANU S1000 Tissue Optimization Kit (BMKMANU, ST03003). Permeabilization, reverse transcription and cDNA synthesis were also performed according to the user guide. Libraries were sequenced on the Illumina NavoSeq. After filtering out low-quality sequences, the remaining sequences were aligned to the chicken reference genome using BSTMatrix v2.3 with default parameters, and an expression profile matrix was generated.

To locate the cell states in the transcriptomics slides, the combined scRNA-seq datasets were integrated with spatial transcriptomics data using conditional autoregressive-based deconvolution (CARD)⁹⁷. This method allowed the transfer of cell-type annotations from scRNA-seq to spatial transcriptomics. CARD was also used to infer correlations in cell-type proportions across spatial locations between pairs of cell types. After annotating cell type diversity within each spot, significant ligand-receptor pairs between neighboring spots were visualized using ggplot2.

Pseudotime trajectory analysis

The cell lineage trajectory was inferred via Monocle2 or Slingshot^98,99. The trajectory of chicken FDCs was constructed by Monocle2. DDRtree was used to visualize in two-dimensional space. Slingshot was used to define computationally imputed pseudotime trajectories of γδ T cells. Double negative T cells were considered as the root state when calculating the trajectories and the pseudotime.

Identifying cell type-associated TFs

To uncover cell type-specific regulation, GENIE3¹⁰⁰ and SCENIC¹⁰¹ were applied to infer the gene regulatory network and calculate the TF activity scores. The regulons specifically active in a selected cluster (compared to other clusters) were identified using the “calcRSS” and “plotRSS” functions in the AUCell R package. After screening the target modules, the gene interaction network was visualized using cytoscape¹⁰² and Gephi¹⁰³.

Cell-cell interaction analysis

The R package CellChat¹⁰⁴ was used to identify and visualize cellular cross-talk between different cell types. The over-expressed ligands or receptors were identified by “identifyOverExpressedGenes” and “identifyOverExpressedInteractions” functions. Then, we used the “computeCommunProb” and “filterCommunication” functions to compute communication probability and infer cellular communication networks. Cell-cell communication at the signaling pathway level between cell types was inferred using the “computeCommunProbPathway” and “aggregateNet” functions.

Spatial co-localization of receptor-ligand pairs

After annotating cell type diversity within each spot, ligand-receptor hotspots were defined as spatial regions containing high numbers of interacting cells and high ligand-receptor co-expression. The expression of receptor-ligand pairs was visualized by Seurat.

Similarity of cell types across amniotes

For evolutionary changes in gene expression, we binarized genes as either expressed or not based on each cell type's average expression profile¹⁹. A gene was considered expressed if it was present in more than 25% of cells within a specific cell type and its expression level was above the third quartile for that cell type. Unsupervised MetaNeighbor analysis was used to systematically assess the transcriptional similarity between cell types across species, with mean AUROC (Area Under the Receiver Operating Characteristic) scores quantifying the similarity of cell-type pairs ³⁴. Before conducting the MetaNeighbor analysis, we quantified expression levels as counts per million (CPM) for each cell in the atlases and then scaled these values using the natural logarithm transformation (ln(CPM + 1)). To address data sparsity in low-coverage sequencing datasets, we employed a pseudo-cell approach to aggregate data from multiple cells of the same cell type. Specifically, we selected 30 cells from each cell type within each species to construct pseudo-cells. The dendrogram of immune cell types was then generated by hierarchical clustering using the Ward.D2 method, with distances based on mean AUROC scores.

Cross-species comparisons of transcriptomics data of immune cells

For cross-species comparisons, we used previously published comprehensive immune cell datasets from humans¹⁰⁵. We then integrated the immune cells of three species using Seurat’s SCTransform workflow. Each cell class (Mononuclear phagocytic cells, B cells, and T cells and innate lymphocytes) was integrated separately. Only one-to-one ortholog genes between these three species, identified using EggNOG mapper and Ensembl, were retained. Integration features were selected by the “SelectIntegrationFeature” function. Integration anchors were identified based on the first 40 canonical components (“FindIntegrationAnchors” function, reduction= “CCA” and normalization.method = “SCT”)¹⁰⁶. The first 18 principal components of the integrated data were used to visualize the data by a UMAP embedding and to construct a neighbourhood graph. For integrated T cells and innate lymphocyte cells, we obtained mean expression profiles for each cluster. The cluster distance matrix was computed as Spearman correlation and used for hierarchical clustering with the Ward.D2 method to generate dendrograms. Dendrograms were coloured according to the proportion of cells from each species within the integrated cluster.

Chicken and human gene alignment

We calculated mean expression profiles for each cell type (scaled as ln(CPM + 1)) and pairwise pearson correlation coefficients using the “cor” function¹⁹. Species-enriched gene expression was defined as genes enriched 2-fold in either direction (chicken > human or human > chicken) with a p-value less than 0.05 (calculated by “MAST”).

Weighted Gene Co-expression Network Analysis

Weighted gene co-expression network analysis (WGCNA) was performed with functions in the WGCNA R package¹⁰⁷. To attenuate the effects of noise and outliers, we aggregated data from 11 to 50 cells within the same cluster to create pseudo-cells for each cell type. The top 2,500 variable genes were used for analysis. The adjacency matrix was constructed by setting the soft power parameter to 10 using “pickSoftThreshold” function. From this adjacency matrix, a topological overlap matrix (TOM) was calculated using the “TOMsimilarityFromExpr” function and the TOM dissimilarity measure (1- TOM) was then used to cluster genes. Modules were identified using the dynamic tree-cutting algorithm with the “cutreeDynamic” function, and module eigengenes were defined using the “moduleEigengene” function.

Gene expression pattern

We categorized expression patterns into four types based on a previously established method ³⁵. Briefly, both type 0 and type 1 involve comparisons of gene expression changes within the same cell type between chickens and humans. Type 0 represents genes expressed in both species, while Type 1 indicates a simple gain or loss of expression between the species. Type 2 changes involve the gain or loss of expression in additional cell types during evolution. Type 3 refers to a change in expression from one cell type to another.

Enrichment analysis for GO and KEGG

Gene Ontology (GO) and KEGG pathways analyses were performed using David (https://david.ncifcrf.gov/summary.jsp) and Metascape (https://metascape.org). GO terms and KEGG pathways with p < 0.05 were considered significantly enriched.

Data availability

The datasets analyzed in this study are available from the Gene Expression Omnibus (GEO) repository under the following accession numbers: PRJNA1125639 and PRJNA1128021.

Human scRNA-seq datasets were collected from the published articles and database (GSE134355; E-MTAB-11536; Gut Cell Atlas: https://www.gutcellatlas.org/).

Competing interests

The authors declare no competing interests.

Author contributions

Y.J., X.F., F.W, and L.F. conceived and designed the project. Y.J., X.F. and L.F., supervised the work. F.W., J.R., Y.Z., H.L. and W.H. collected tissue samples and performed data analyses. Y.J., M.Q., T.S., H.S., H.T., H.W. and X.G. helped with bioinformatic analyses. Y.J., X.H., X.F., and Z.H. provided funding. J.M., Z.Y., L.F., Y.J., and H.Z. provided valuable comments. F.W., J.R. and Y.Z. wrote the manuscript. Y.J., L.F., J.S., H.Z., M.F., S.E.D., L.L. and Y.L. revised the manuscript. All authors contributed to the work. All authors read and approved the manuscript for submission.

Acknowledgments

This work was supported by the National Key R&D Program of China (2022YFF1000100 to Y.J.), China Agriculture Research System of MOF and MARA (CARS-41-Z08 to H.Z.), Hunan Poultry Industry Technology System, the National Natural Science Foundation of China (32302733 to F.W.), China Postdoctoral Science Foundation (313974 to F.W.), Postdoctoral Research Project of Shaanxi Province (2023BSHYDZZ86 to F.W.). We thank the High-Performance Computing platform of Northwest A&F University.

Benton, M. J. & Donoghue, P. C. J., Paleontological evidence to date the tree of life. MOL BIOL EVOL24 26 (2007).
Brown, W. R. A., Hubbard, S. J., Tickle, C. & Wilson, S. A., The chicken as a model for large-scale analysis of vertebrate gene function. NAT REV GENET4 87 (2003).
Morgan, B. A., Izpisúa-Belmonte, J. C., Duboule, D. & Tabin, C. J., Targeted misexpression of Hox-4.6 in the avian limb bud causes apparent homeotic transformations. NATURE358 236 (1992).
Olsen, B. et al., Global patterns of influenza a virus in wild birds. SCIENCE312 384 (2006).
Smith, J. C. et al., Prevalence and molecular characterization of Salmonella isolated from wild birds in fresh produce environments. FRONT MICROBIOL14 1272916 (2023).
Morera, D. et al., RNA-Seq reveals an integrated immune response in nucleated erythrocytes. PLOS ONE6 e26998 (2011).
Oláh, I. & Nagy, N., Retrospection to discovery of bursal function and recognition of avian dendritic cells; past and present. DEV COMP IMMUNOL41 310 (2013).
Cihak, J. et al., In vivo depletion of chicken T-cell subsets. SCAND J IMMUNOL38 123 (1993).
Han, X. et al., Construction of a human cell landscape at single-cell level. NATURE581 303 (2020).
Han, L. et al., Cell transcriptomic atlas of the non-human primate Macaca fascicularis. NATURE604 723 (2022).
Wang, F. et al., Endothelial cell heterogeneity and microglia regulons revealed by a pig cell landscape at single-cell level. NAT COMMUN13 3620 (2022).
Maxwell, M., Söderlund, R., Härtle, S. & Wattrang, E., Single-cell RNA-seq mapping of chicken peripheral blood leukocytes. BMC GENOMICS25 124 (2024).
Cheng, J. et al., B Lymphocyte Development in the Bursa of Fabricius of Young Broilers is Influenced by the Gut Microbiota. MICROBIOL SPECTR11 e0479922 (2023).
Dai, M. et al., Dissection of key factors correlating with H5N1 avian influenza virus driven inflammatory lung injury of chicken identified by single-cell analysis. PLOS PATHOG19 e1011685 (2023).
Yamagata, M., Yan, W. & Sanes, J. R., A cell atlas of the chick retina based on single-cell transcriptomics. ELIFE10 (2021).
Mantri, M. et al., Spatiotemporal single-cell RNA sequencing of developing chicken hearts identifies interplay between cellular differentiation and morphogenesis. NAT COMMUN12 1771 (2021).
HAMBURGER, V. & HAMILTON, H. L., A series of normal stages in the development of the chick embryo. J MORPHOL88 49 (1951).
Vidarsson, H. et al., The forkhead transcription factor Foxi1 is a master regulator of vacuolar H-ATPase proton pump subunits in the inner ear, kidney and epididymis. PLOS ONE4 e4471 (2009).
Travaglini, K. J. et al., A molecular cell atlas of the human lung from single-cell RNA sequencing. NATURE587 619 (2020).
Nombela, I., Ortega-Villaizan, M. D. M. & Spindler, K. R., Nucleated red blood cells: Immune cell mediators of the antiviral response. PLOS PATHOG14 e1006910 (2018).
Burgoyne, R. D., Mechanisms of catecholamine secretion from adrenal chromaffin cells. J PHYSIOL PHARMACOL46 273 (1995).
Voets, T. et al., Intracellular calcium dependence of large dense-core vesicle exocytosis in the absence of synaptotagmin I. P NATL ACAD SCI USA98 11680 (2001).
Hepp, R. et al., Adrenal gland SNAP-25 expression is altered in thyroid hormone receptor knock-out mice. NEUROREPORT12 1427 (2001).
Malaiyandi, L. M. et al., M(3)-subtype muscarinic receptor activation stimulates intracellular calcium oscillations and aldosterone production in human adrenocortical HAC15 cells. MOL CELL ENDOCRINOL478 1 (2018).
Zhang, Y. et al., Growth hormone (GH) receptor (GHR)-specific inhibition of GH-Induced signaling by soluble IGF-1 receptor (sol IGF-1R). MOL CELL ENDOCRINOL492 110445 (2019).
Dong, J. et al., ACACA reduces lipid accumulation through dual regulation of lipid metabolism and mitochondrial function via AMPK- PPARα- CPT1A axis. J TRANSL MED22 196 (2024).
Dai, W. et al., FASN deficiency induces a cytosol-to-mitochondria citrate flux to mitigate detachment-induced oxidative stress. CELL REP42 112971 (2023).
Mullins, L. J. et al., Mineralocorticoid Excess or Glucocorticoid Insufficiency: Renal and Metabolic Phenotypes in a Rat Hsd11b2 Knockout Model. HYPERTENSION66 667 (2015).
Han, F. et al., HSD3B1 variant and androgen-deprivation therapy outcome in prostate cancer. CANCER CHEMOTH PHARM87 103 (2021).
Lee, S. C., Baranowski, E. S., Sakremath, R., Saraff, V. & Mohamed, Z., Hypoglycaemia in adrenal insufficiency. FRONT ENDOCRINOL14 1198519 (2023).
Lerario, A. M., Finco, I., LaPensee, C. & Hammer, G. D., Molecular Mechanisms of Stem/Progenitor Cell Maintenance in the Adrenal Cortex. FRONT ENDOCRINOL8 52 (2017).
de Mendonca, P. O. R., Costa, I. C. & Lotfi, C. F. P., The involvement of Nek2 and Notch in the proliferation of rat adrenal cortex triggered by POMC-derived peptides. PLOS ONE9 e108657 (2014).
Iwahashi, N. et al., Single-cell and spatial transcriptomics analysis of human adrenal aging. MOL METAB84 101954 (2024).
Fischer, S., Crow, M., Harris, B. D. & Gillis, J., Scaling up reproducible research for single-cell transcriptomics using MetaNeighbor. NAT PROTOC16 4031 (2021).
Travaglini, K. J. et al., A molecular cell atlas of the human lung from single-cell RNA sequencing. NATURE587 619 (2020).
Gilliland, H. N., Beckman, O. K. & Olive, A. J., A Genome-Wide Screen in Macrophages Defines Host Genes Regulating the Uptake of Mycobacterium abscessus. MSPHERE8 e0066322 (2023).
Nagy, N., Magyar, A., Dávid, C., Gumati, M. K. & Oláh, I., Development of the follicle-associated epithelium and the secretory dendritic cell in the bursa of fabricius of the guinea fowl (Numida meleagris) studied by novel monoclonal antibodies. The Anatomical record262 279 (2001).
Nagy, N., Bódi, I. & Oláh, I., Avian dendritic cells: Phenotype and ontogeny in lymphoid organs. DEV COMP IMMUNOL58 47 (2016).
Malaguarnera, L., Marsullo, A., Zorena, K., Musumeci, G. & Di Rosa, M., Vitamin D(3) regulates LAMP3 expression in monocyte derived dendritic cells. CELL IMMUNOL311 13 (2017).
Hontelez, S. et al., Dendritic cell-specific transcript: dendritic cell marker and regulator of TLR-induced cytokine production. J IMMUNOL189 138 (2012).
Haldar, M. et al., Heme-mediated SPI-C induction promotes monocyte differentiation into iron-recycling macrophages. CELL156 1223 (2014).
Kelly, L. M., Englmeier, U., Lafon, I., Sieweke, M. H. & Graf, T., MafB is an inducer of monocytic differentiation. The EMBO journal19 1987 (2000).
Ansel, K. M. et al., A chemokine-driven positive feedback loop organizes lymphoid follicles. NATURE406 309 (2000).
Badr, G. et al., BAFF enhances chemotaxis of primary human B cells: a particular synergy between BAFF and CXCL13 on memory B cells. BLOOD111 2744 (2008).
Ratcliffe, M. J. H. & Härtle, S., in Avian Immunology (Third Edition), edited by Bernd Kaspers, Karel A. Schat, Thomas W. Göbel, and Lonneke Vervelde (Academic Press, 2022), pp. 71.
Mackay, F. & Browning, J. L., BAFF: a fundamental survival factor for B cells. NAT REV IMMUNOL2 465 (2002).
Suo, C. et al., Mapping the developing human immune system across organs. SCIENCE376 eabo0510 (2022).
Heesters, B. A., Myers, R. C. & Carroll, M. C., Follicular dendritic cells: dynamic antigen libraries. NAT REV IMMUNOL14 495 (2014).
Rodda, L. B. et al., Single-Cell RNA Sequencing of Lymph Node Stromal Cells Reveals Niche-Associated Heterogeneity. IMMUNITY48 1014 (2018).
Matz, H. & Dooley, H., 450 million years in the making: mapping the evolutionary foundations of germinal centers. FRONT IMMUNOL14 1245704 (2023).
Matz, H. et al., Organized B cell sites in cartilaginous fishes reveal the evolutionary foundation of germinal centers. CELL REP42 112664 (2023).
Wittner, J. & Schuh, W., Krüppel-like factor 2: a central regulator of B cell differentiation and plasma cell homing. FRONT IMMUNOL14 1172641 (2023).
Clottu, A. S. et al., EBI2 Expression and Function: Robust in Memory Lymphocytes and Increased by Natalizumab in Multiple Sclerosis. CELL REP18 213 (2017).
Kurosaki, T., Kometani, K. & Ise, W., Memory B cells. NAT REV IMMUNOL15 149 (2015).
Laidlaw, B. J., Duan, L., Xu, Y., Vazquez, S. E. & Cyster, J. G., The transcription factor Hhex cooperates with the corepressor Tle3 to promote memory B cell development. NAT IMMUNOL21 1082 (2020).
Kumar, S. & Hedges, S. B., A molecular timescale for vertebrate evolution. NATURE392 917 (1998).
Muramatsu, M. et al., Class switch recombination and hypermutation require activation-induced cytidine deaminase (AID), a potential RNA editing enzyme. CELL102 553 (2000).
GREY, H. M., Phylogeny of the immune response. studies on some physical chemical and serologic characteristics of antibody produced in the turtle. J IMMUNOL91 819 (1963).
Marchalonis, J. J., Ealey, E. H. & Diener, E., Immune response of the Tuatara, Sphenodon punctatum. The Australian journal of experimental biology and medical science47 367 (1969).
Liao, Y. et al., Cell landscape of larval and adult Xenopus laevis at single-cell resolution. NAT COMMUN13 4306 (2022).
Lin, L., Gerth, A. J. & Peng, S. L., Active inhibition of plasma cell development in resting B cells by microphthalmia-associated transcription factor. The Journal of experimental medicine200 115 (2004).
Yamanashi, Y. et al., Role of the rasGAP-associated docking protein p62(dok) in negative regulation of B cell receptor-mediated signaling. GENE DEV14 11 (2000).
Buchta, C. M. & Bishop, G. A., TRAF5 negatively regulates TLR signaling in B lymphocytes. J IMMUNOL192 145 (2014).
Porstner, M. et al., miR-148a promotes plasma cell differentiation and targets the germinal center transcription factors Mitf and Bach2. EUR J IMMUNOL45 1206 (2015).
Capasso, M. et al., HVCN1 modulates BCR signal strength via regulation of BCR-dependent generation of reactive oxygen species. NAT IMMUNOL11 265 (2010).
Arendt, D. et al., The origin and evolution of cell types. NAT REV GENET17 744 (2016).
Maxwell, M., Söderlund, R., Härtle, S. & Wattrang, E., Single-cell RNA-seq mapping of chicken peripheral blood leukocytes. BMC GENOMICS25 124 (2024).
Hamey, F. K. et al., Single-cell molecular profiling provides a high-resolution map of basophil and mast cell development. ALLERGY76 1731 (2021).
Ribot, J. C., Lopes, N. & Silva-Santos, B., γδ T cells in tissue physiology and surveillance. NAT REV IMMUNOL21 221 (2021).
Malhotra, N. et al., A network of high-mobility group box transcription factors programs innate interleukin-17 production. IMMUNITY38 681 (2013).
Zuberbuehler, M. K. et al., The transcription factor c-Maf is essential for the commitment of IL-17-producing γδ T cells. NAT IMMUNOL20 73 (2019).
Bonneville, M., O'Brien, R. L. & Born, W. K., Gammadelta T cell effector functions: a blend of innate programming and acquired plasticity. NAT REV IMMUNOL10 467 (2010).
Shanebeck, K. D. et al., Regulation of murine B cell growth and differentiation by CD30 ligand. EUR J IMMUNOL25 2147 (1995).
Mitsdoerffer, M. et al., Proinflammatory T helper type 17 cells are effective B-cell helpers. P NATL ACAD SCI USA107 14292 (2010).
Grogan, J. L. & Ouyang, W., A role for Th17 cells in the regulation of tertiary lymphoid follicles. EUR J IMMUNOL42 2255 (2012).
Qu, G. et al., Comparing Mouse and Human Tissue-Resident γδ T Cells. FRONT IMMUNOL13 891687 (2022).
Tosches, M. A. et al., Evolution of pallium, hippocampus, and cortical cell types revealed by single-cell transcriptomics in reptiles. SCIENCE360 881 (2018).
Hain, D. et al., Molecular diversity and evolution of neuron types in the amniote brain. SCIENCE377 eabp8202 (2022).
Murat, F. et al., The molecular evolution of spermatogenesis across mammals. NATURE613 308 (2023).
Carroll, S. B., Evo-devo and an expanding evolutionary synthesis: a genetic theory of morphological evolution. CELL134 25 (2008).
Jiao, A. et al., Single-cell sequencing reveals the evolution of immune molecules across multiple vertebrate species. J ADV RES55 73 (2024).
Krautler, N. J. et al., Follicular dendritic cells emerge from ubiquitous perivascular precursors. CELL150 194 (2012).
Aguzzi, A. & Krautler, N. J., Characterizing follicular dendritic cells: A progress report. EUR J IMMUNOL40 2134 (2010).
Tew, J. G., Wu, J., Fakher, M., Szakal, A. K. & Qin, D., Follicular dendritic cells: beyond the necessity of T-cell help. TRENDS IMMUNOL22 361 (2001).
Vainio, O., Viljanen, M. K. & Toivanen, A., Early ontogeny of germinal center formation in the chicken. DEV COMP IMMUNOL2 493 (1978).
Schussek, S. et al., The CTA1-DD adjuvant strongly potentiates follicular dendritic cell function and germinal center formation, which results in improved neonatal immunization. MUCOSAL IMMUNOL13 545 (2020).
Hoshi, H. et al., Patterns of age-dependent changes in the numbers of lymph follicles and germinal centres in somatic and mesenteric lymph nodes in growing C57Bl/6 mice. J ANAT198 189 (2001).
Laursen, A. M. S. et al., Characterizaton of gamma delta T cells in Marek's disease virus (Gallid herpesvirus 2) infection of chickens. VIROLOGY522 56 (2018).
Berndt, A., Pieper, J. & Methner, U., Circulating gamma delta T cells in response to Salmonella enterica serovar enteritidis exposure in chickens. INFECT IMMUN74 3967 (2006).
Hein, W. R. & Mackay, C. R., Prominence of gamma delta T cells in the ruminant immune system. Immunology today12 30 (1991).
Holderness, J., Hedges, J. F., Ramstead, A. & Jutila, M. A., Comparative biology of γδ T cell function in humans, mice, and domestic animals. ANNU REV ANIM BIOSCI1 99 (2013).
Zheng, G. X. Y. et al., Massively parallel digital transcriptional profiling of single cells. NAT COMMUN8 14049 (2017).
Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R., Integrating single-cell transcriptomic data across different conditions, technologies, and species. NAT BIOTECHNOL36 411 (2018).
Korsunsky, I. et al., Fast, sensitive and accurate integration of single-cell data with Harmony. NAT METHODS16 1289 (2019).
Chen, A. et al., Spatiotemporal transcriptomic atlas of mouse organogenesis using DNA nanoball-patterned arrays. CELL185 1777 (2022).
Liu, H. et al., Integrative molecular and spatial analysis reveals evolutionary dynamics and tumor-immune interplay of in situ and invasive acral melanoma. CANCER CELL42 1067 (2024).
Ma, Y. & Zhou, X., Spatially informed cell-type deconvolution for spatial transcriptomics. NAT BIOTECHNOL40 1349 (2022).
Trapnell, C. et al., The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. NAT BIOTECHNOL32 381 (2014).
Street, K. et al., Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC GENOMICS19 477 (2018).
Huynh-Thu, V. A., Irrthum, A., Wehenkel, L. & Geurts, P., Inferring regulatory networks from expression data using tree-based methods. PLOS ONE5 (2010).
Aibar, S. et al., SCENIC: single-cell regulatory network inference and clustering. NAT METHODS14 1083 (2017).
Shannon, P. et al., Cytoscape: a software environment for integrated models of biomolecular interaction networks. GENOME RES13 2498 (2003).
Bastian, M., Heymann, S. & Jacomy, M., Gephi: An Open Source Software for Exploring and Manipulating Networks. Proceedings of the International AAAI Conference on Web and Social Media3 361 (2009).
Jin, S. et al., Inference and analysis of cell-cell communication using CellChat. NAT COMMUN12 1088 (2021).
Domínguez Conde, C. et al., Cross-tissue immune cell analysis reveals tissue-specific features in humans. SCIENCE376 eabl5197 (2022).
Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R., Integrating single-cell transcriptomic data across different conditions, technologies, and species. NAT BIOTECHNOL36 411 (2018).
Langfelder, P. & Horvath, S., WGCNA: an R package for weighted correlation network analysis. BMC BIOINFORMATICS9 559 (2008).

There is NO Competing Interest.

Download PDF

Version 1

posted

You are reading this latest preprint version

Comparison of single-cell landscapes among three amniotes reveals conservation and innovation in avian immune systems

Status:

Version 1

Abstract

Figures

Introduction

Results

Construction of single-cell and spatial transcriptome atlas for three amniotes

Comparison of cell landscapes among chickens, turtles and humans

Cellular divergence driven by evolutionary lability in gene expression localization

The evolution of GC B cells and affinity maturation of antibodies in amniotes

Divergent evolution of the γδ T lineage in amniotes

Avian γδ T cell function and ontogeny

Discussion

Methods

Ethics statement

Collection of chicken, duck and turtle tissues

Single-nucleus/cell suspension preparation

scRNA-seq and snRNA-seq data processing

Spatial transcriptomics sequencing and data analysis

Pseudotime trajectory analysis

Identifying cell type-associated TFs

Cell-cell interaction analysis

Spatial co-localization of receptor-ligand pairs

Similarity of cell types across amniotes

Cross-species comparisons of transcriptomics data of immune cells

Chicken and human gene alignment

Weighted Gene Co-expression Network Analysis

Gene expression pattern

Enrichment analysis for GO and KEGG

Declarations

Data availability

Competing interests

Author contributions

Acknowledgments

References

Additional Declarations

Supplementary Files

Status:

Version 1