Tree shrew immune cell atlas identifies NR1H3⁺ tissue macrophages with conserved anti-inflammatory function

Domain: neuroscience
Published: 2026-04-09
Editorial outcome: major_revision

Abstract
Rodents are widely used in immunology but do not always recapitulate human immune functions. The tree shrew (Tupaia belangeri) is phylogenetically closer to primates than rodents and may help bridge this gap, yet its immune system has not been comprehensively characterised at single-cell resolution. Here, we present a single-cell transcriptomic atlas of the tree shrew immune system, profiling 39 cell types across 12 tissues. We uncover human-like tonsillar structures and two transcriptionally distinct splenic macrophage subsets: an NR1H3+ immunoregulatory subset and an IRF8+ pro-inflammatory subset. During acute Epstein-Barr virus infection, NR1H3+ macrophages undergo pronounced inflammatory reprogramming and emerge as a major intercellular signalling hub. Cross-species integration of single-cell datasets from seven vertebrate species indicates that NR1H3+ macrophage gene programs are evolutionarily conserved. Functional analyses further support a shared anti-inflammatory programme mediated by NR1H3-dependent suppression of the non-canonical NF-κB pathway in both tree shrew and human macrophages. Together, these findings provide a reference resource for tree shrew immunology and support the tree shrew as a complementary model alongside rodent and primate systems.

Introduction
Animal models play indispensable roles in understanding human diseases and developing effective therapeutic interventions. Rodents are widely used because they are accessible and supported by well-established experimental tools, but they do not always recapitulate human physiology and pathology, raising concerns about translational studies1–3. For example, mice and rats are not naturally susceptible to hepatitis viruses4, herpesviruses5 or SARS-CoV-26, partly because they lack key entry receptors, which limits their utility for studying these infections.

Nonhuman primates (NHPs) can provide more translatable insights. For example, smaller New World NHPs such as cottontop tamarins and common marmosets have served as experimental systems for Epstein–Barr virus (EBV) positive lymphoma models7. However, ethical, cost and availability constraints limit NHPs research, underscoring the need for alternative small-animal models that more closely mirror human biology8. One promising candidate is the tree shrew, a small mammal that combines the practical advantages of rodents (compact size, rapid reproduction, short lifespan, and cost-effectiveness) with a closer evolutionary affinity to primates9,10. Tree shrews share many physiological and molecular features with primates, including immune system components11,12, and are susceptible to several human pathogens that do not infect mice, such as hepatitis B/C/D viruses, EBV and Zika virus13–16. Genomic analyses further support their close phylogenetic relationship with primates17,18. Despite these advantages, the paucity of species-specific research tools (e.g., well-characterized antibodies) has hampered broader application of tree shrews, particularly in immunological studies that rely on defined cell markers19.

Advances in single-cell genomics now offer a way to surmount this limitation. High-throughput single-cell RNA sequencing (scRNA-seq) technology enables comprehensive characterization of cell populations without the need for species-specific reagents20. This approach provides an opportunity to characterize nontraditional model organisms such as tree shrew at molecular resolution. Indeed, recent cross-species single-cell studies have uncovered immune cell subsets and gene regulatory programs, illustrating the value of this approach21–25.

Here we map the immune landscape of the tree shrew using single-cell and single-nucleus RNA sequencing across 12 tissues. We identify 39 distinct cell types and reveal well-developed tonsillar structures analogous to human tonsils, which are absent in rodents. Moreover, we define two splenic macrophage states (NR1H3+ and IRF8+) with distinct inflammatory programs that shift during acute EBV infection. Cross-species integration and functional assays show that NR1H3+ macrophages and NR1H3-dependent suppression of non-canonical NF-κB signalling are conserved between tree shrew and human. Together, these findings provide a reference for tree shrew immunology and support its use as a complementary small-animal model for translational and comparative studies.

Results
Construction of a multitissue single-cell atlas reveals the cellular diversity of tree shrew

To systematically characterize the cellular heterogeneity of the tree shrew, we constructed a high-resolution single-cell atlas using scRNA-seq (Fig. 1A). We generated 24 sc/snRNA‑seq libraries from 12 tissues from three adult tree shrews. Immune‑related tissues were sampled across animals, whereas non‑immune parenchymal tissues and the pharyngeal tonsil were profiled once because of scope and tissue availability (Fig. 1B and Supplementary Data 1). Most tissues were analysed via scRNA-seq, but single-nucleus RNA sequencing (snRNA-seq) was used for the heart and liver to mitigate dissociation-associated bias for large or fragile cells26. A total of 221,113 cells/nuclei were recovered and passed initial quality control, and 186,973 high-quality cells/nuclei were retained after downstream filtering and doublet removal. Cell numbers ranged from 2,765 in pharyngeal tonsil to 30,645 in lymph tissue (Fig. 1B). We integrated all of the scRNA-seq/snRNA-seq data using the reciprocal principal component analysis (RPCA) method to correct for potential batch effects27. Unsupervised clustering and uniform manifold approximation and projection (UMAP) clearly delineated major cell types, including endothelial cells, fibroblasts, haematopoietic stem cells and multipotent progenitor cells (HSC/MPP), myeloid cells, parenchymal cells, platelets, B and T cells (Supplementary Fig. 1). Given cost considerations, we prioritised biological replication for immune tissues; most non-immune parenchymal tissues were profiled once. Nevertheless, cells from the same tissue (or functionally related immune tissues such as lymph nodes, spleen and palatine tonsils) showed consistent distributions across animals (Supplementary Fig. 2), supporting the robustness of our atlas. Moreover, profiling multiple tissues from the same individual allowed us to control for inter-individual variables such as age, sex, diet and environment.

Further high-resolution clustering based on the expression of marker genes identified 39 distinct cell clusters within eight major cell classes in the global UMAP of all tissues (Fig. 1C and Supplementary Fig. 3). On average, each cell expressed 1,985 genes and 6,641 unique molecular identifiers (UMIs). The number of cells varied across the 39 cell types, ranging from 128 proliferating T cells to 17,776 double-positive T cells (DP T) (Supplementary Fig. 4). As expected, parenchymal subclusters were largely tissue-restricted, whereas multiple immune cell types were shared across tissues (Supplementary Fig. 5). Tissue-preference analysis further indicated that each tissue has a distinct immune neighbourhood; for example, although spleen and lymph nodes were both enriched for B cells, their myeloid compositions differed (Fig. 1D).

To further validate atlas annotations and explore underlying regulatory programs, we performed a SCENIC analysis to reconstruct transcription factor regulons across cells/ tissues (Supplementary Fig. 6 and 7). Each module was driven by transcription factors whose activity was enriched in particular cell populations or tissues. These regulon modules showed expression patterns consistent with known cell-type distribution, and their target genes were enriched for functions characteristic of those cells or tissues. For instance, modules active in metabolically focused parenchymal cells were enriched for metabolic pathways, whereas those prominent in proliferative progenitors were enriched in cell-cycle processes. The analysis also highlighted key lineage-defining TFs (e.g. HNF1B in kidney cells28 and FOXM1/E2F-family in haematopoietic progenitors29), supporting the accuracy of our cell-type annotations (Supplementary Fig. 8 and 9). In an independent bone-marrow analysis, regulon-based clustering recapitulated the transcriptome-defined cell populations (Supplementary Fig. 10), further confirming the robustness of our cell identification. Overall, this integrative transcription factor network analysis confirmed that each identified cell type possesses a distinct and biologically coherent regulatory program.

Anatomical and functional validation of tonsillar immunity in tree shrew

To further investigate intratissue cellular heterogeneity, we focused on the tonsils, which lie at the intersection of the digestive respiratory tracts30. Unlike rodents, which rely on nasal-associated lymphoid tissue as the primary upper airway immune organ, tree shrews possess human-like tonsillar structures, offering a useful model for mucosal immunity31,32. Coronal head sections revealed two distinct lymphoid structures: palatine tonsils in the oropharyngeal mucosa and pharyngeal tonsils in the pharyngeal region (Fig. 2A). Single-cell sequencing of the palatine tonsils identified 17 distinct cell populations span immune and epithelial/stromal compartments, including B cells (CD79A/MS4A1; plasma cells MZB1/JCHAIN), CD4⁺ T‑cell states (naïve LEF1/SELL, memory S100A4/TNFRSF18, Treg FOXP3/LAG3, exhausted CTLA4/TOX), myeloid populations (macrophages CSF1R/CD68, dendritic cells FCER1A/ITGAX), and epithelial/secretory cells (EPCAM/KRT14/KRT16) (Fig. 2B, C). Pharyngeal tonsil cells were analysed separately (Supplementary Fig. 11). These markers are consistent with those reported for the corresponding human cell types33. For example, LEF is required for T-cell development, and FOXP3 is indispensable for regulatory T-cell suppressive function34,35. KEGG pathway analysis further revealed cell type-specific signatures. For example, the B- and T-cell receptor signalling pathways are enriched in B- and T-cell subtypes, respectively. Pathways such as the cytokine‒cytokine receptor interaction pathway, JAK‒STAT signalling pathway, and PI3K‒Akt signalling pathway were enriched in NK cells. Plasma cells were enriched for protein processing in the endoplasmic reticulum, consistent with their role in immunoglobulin production36 (Fig. 2D). We subsequently validated representative markers in palatine tonsils by immunohistochemistry. CD3D (T-lineage) was prominent in peripheral regions of germinal centres, whereas CD79A (B-lineage) showed strong staining within germinal centres (Fig. 2E). MKI67, which marks proliferating cells across active phases of the cell cycle37, was detected in basal epithelial cells and germinal-centre cells. BCL6 and CR2 were enriched in germinal-centre cells, consistent with findings from human tonsil scRNA-seq33. Finally, KRT14 marked epithelial cells.

In tree shrew pharyngeal tonsils, ciliated cells expressed TPPP3, NEK10, and EPCAM, whereas glandular cells expressed AQP5 and FAM3D. Goblet cells expressed KRT15 and VMO1. Immune cells, including B cells, T cells and macrophages, were also present (Supplementary Fig. 11A-D). KEGG enrichment analysis showed enrichment of B- and T-cell receptor signalling pathways in B and T cells, respectively. Ciliated cells were enriched for thiamine metabolism, consistent with high energy demand, whereas glandular cells were enriched for salivary secretion pathways38 (Supplementary Fig. 11E). Immunohistochemical staining confirmed canonical markers (CD3D, CD79A, MKI67, BCL6, CR2 and AQP5) in pharyngeal tonsil sections (Supplementary Fig. 11F). Although scRNA-seq recovered relatively few B cells and plasma cells and did not capture mature/activated B-cell subsets, immunohistochemistry revealed germinal centres, indicating that these populations are present in vivo and were likely under-sampled because of limited tissue yield.

Tissue-specific functional diversity of tree shrew macrophages

Following comprehensive characterization of baseline cellular diversity, we next investigated the tissue-specific adaptation of macrophages in tree shrew. Although tissue-resident macrophages are recognized as pivotal regulators of tissue homeostasis and inflammatory responses in mammals39,40, their functional specialization and transcriptional landscape remain unexplored in this emerging model organism. We extracted 10,466 high-quality cells previously annotated as macrophages from all 12 tissues. Expression of the conserved macrophage marker CSF1R confirmed their identity (Fig. 3A). Immunofluorescence staining in spleen tissue further validated CSF1R expression, with CSF1R+ macrophages anatomically localised at the periphery of BCL6+ germinal centres (Fig. 3B), consistent with expected localisation. Macrophages from stomach (n = 44), pharyngeal tonsil (n = 92) and kidney (n = 0) were excluded because of insufficient cell numbers for robust subpopulation analysis. After batch correction, unsupervised clustering revealed nine transcriptionally distinct macrophage subpopulations (Fig. 3C, Supplementary Data 2). Tissue preference analysis using observed-to-expected cell ratios (Ro/e)41 demonstrated striking tissue specificity, with each subcluster showing strong enrichment in particular tissues (Ro/e >5) (Fig. 3D). Functional annotation of subcluster-defining genes (top 50 markers) by Gene Ontology (GO) analysis provided evidence for tissue microenvironment-driven specialization (Fig. 3E). For example, hepatic macrophages (C2) were enriched in metabolic regulation pathways, including the urea cycle, nitrogen cycle metabolic process, and L-arginine transmembrane transport, consistent with liver detoxification functions. Cardiac macrophages (C3) showed signatures of cardiac muscle contraction, heart morphogenesis, and actin cytoskeleton organization, suggesting roles in myocardial maintenance. In contrast to macrophages from other tissues, spleen macrophages presented two transcriptionally divergent subclusters (C4, C7) that were strongly spleen-enriched (Ro/e > 5) yet displayed functional divergence beyond typical tissue-driven specialization (Fig. 3C, D). Consistent spatial distributions of the C4 and C7 subclusters across triplicate spleen samples in UMAP plots (Supplementary Fig. 12) suggested genuine heterogeneity rather than batch effects. C4 macrophages were characterized by elevated expression of NR1H3 and PPARG (anti-inflammatory markers), whereas C7 macrophages expressed high levels of FLT3, IRF8, and CD226 (proinflammatory markers) (Fig. 3E). GO analysis of C4-defining genes highlighted pathways related to immunosuppression (e.g., negative regulation of macrophage activation, complement activation, and negative regulation of the type II interferon response), whereas C7 markers were linked to leukocyte activation, cell killing, and positive regulation of type II interferon production (Fig. 3E).

NR1H3⁺ macrophages undergo reprogramming during acute EBV infection

To investigate how splenic macrophage subsets respond to inflammatory stress, we leveraged our tree shrew EBV infection model13 and analysed acute infection (3–7 days post-infection, dpi). We focused on the acute phase to minimisNe confounding effects from monocyte-derived macrophage infiltration that is typically observed during chronic infection42, thereby preserving our analytical focus on tissue-resident macrophage populations. After quality control, we obtained 70,200 high-quality splenic scRNA-seq profiles (control: 17,372; 3 dpi: 28,590; 7 dpi: 24,238) that clustered into 10 populations after batch correction (Fig. 4A, B). UMAP visualization revealed a conserved cellular architecture without infection-induced clusters (Supplementary Fig. 13A), enabling direct comparison of preidentified macrophage subsets. A brief T/NKT analysis is provided in Supplementary Fig. 14. CellChat analysis showed that the number of inferred interactions was conserved across groups, whereas overall communication strength increased during acute infection, consistent with intensified intercellular signalling (Fig. 4C). Subcluster-specific communication strength analysis revealed NR1H3_Macro as a major signalling hub and cytokine source, with sharply increased incoming and outgoing interaction strength at 3 dpi (Fig. 4D). By contrast, IRF8_Macro showed reduced communication strength (Fig. 4D). Pathway analysis indicated activation of IFNG signalling (Supplementary Fig. 13B), with IFNG_NKT cells acting as senders and NR1H3_Macro as receivers via IFNGR1/2 (Supplementary Fig. 13D). NR1H3_Macro also transitioned towards a TNF-secreting phenotype and emerged as a key sender through TNF–TNFRSF1B interactions (Supplementary Fig. 13C). Although the abundance of NR1H3_Macro and IRF8_Macro remained stable across infection stages (Fig. 4E, F), their transcriptional states diverged. UCell scoring using subset-specific markers showed attenuation of anti-inflammatory signatures in NR1H3_Macro, whereas IRF8_Macro scores unchanged (Fig. 4G). Consistently, differential expression analysis within NR1H3_Macro revealed infection-induced downregulation of NR1H3 together with upregulation of proinflammatory mediators (IL1B, TNF, TNFSF9 and TNFSF14) and complement components (C1QA, C1QB and C1QC) (Fig. 4H, I). These results suggest that NR1H3 acts as a rheostat that constrains macrophage inflammatory potential, and its suppression permits inflammatory activation.

Conservation of NR1H3-mediated immunosuppression in human macrophages

Given the prominence of NR1H3⁺ macrophages in tree shrew, we asked whether a similar population exists in humans and exhibits conserved functions. NR1H3+ macrophages have been identified in human spleen scRNA-seq datasets43, and their expression profiles correlated with those of tree shrew (r=0.769, p < 0.01) (Fig. 5A, B, Supplementary Fig. 15), supporting conserved functions. To test function, we generated human NR1H3-overexpressing (NR1H3-OE) and NR1H3-knockdown (NR1H3-KO) THP-1 cell lines using lentivirus transduction system (Fig. 5C). Both NR1H3-OE and NR1H3-KO THP-1 cells differentiated normally, indicating that NR1H3 gain or NR1H3 loss of function did not grossly alter THP-1 macrophage differentiation (Fig. 5D). In functional assays, NR1H3 gain in THP‑1 macrophages reduced secretion of IL‑1β, IL‑6 and TNF‑α, modestly decreased uptake of bacterial particles, and diminished IFN‑γ production and effector polarization of primary CD8⁺ T cells while leaving T‑cell proliferation unchanged (Supplementary Fig. 17A–F). Furthermore, pharmacologic activation of NR1H3 with T0901317 in primary human monocyte‑derived macrophages similarly suppressed LPS‑induced IL‑1β, TNF‑α and IL‑6 release (Supplementary Fig. 17G, H), underscoring the translational relevance of NR1H3‑mediated immunoregulation.

We next performed bulk RNA-seq on NR1H3-OE and wild-type (WT) cells. Consistent with findings in the tree shrew model, TNF (log FoldChange = -1.73, p < 0.01), TNFSF14 (log FoldChange = -0.88, p < 0.01), and C1QA (log FoldChange = -1.37, p = 0.04) were downregulated in NR1H3-OE macrophages (Fig. 5E and Supplementary Data 3). Furthermore, key non-canonical NF-κB pathway-associated genes, including NFKBIZ (log FoldChange = -0.77, p = 0.02), NFKB2 (log FoldChange = -1.29, p < 0.01), and MAP3K14 (log FoldChange = -1.36, p < 0.01), were also repressed (Fig. 5E). Gene set enrichment analysis (GSEA) using the Reactome database showed negative enrichment of REACTOME_TNF_RECEPTOR_SUPERFAMILY_TNFSF_MEMBERS_MEDIATING_NON_CANONICAL_NF_KB_PATHWAY in NR1H3-OE macrophages (NES=-1.49, p = 0.047) (Fig. 5F and Supplementary Data 4). Altered expression of selected genes was validated by qRT–PCR (Fig. 5G). To test whether NR1H3 directly binds to these loci, we performed NR1H3 cleavage under targets and tagmentation (Cut&Tag) analysis (Fig. 5H and Supplementary Data 5). At NR1H3-OE-induced gene loci, NR1H3 binding was detected upstream of TNFSF14, NFKB2 and NFKBIZ (Fig. 5I). We further evaluated pharmacological activation using two NR1H3 agonists (T0901317 and GW3965) 44,45. These agonists increased NR1H3 protein abundance in THP-1 macrophages (Fig. 5J). Consistent with NR1H3-mediated repression, agonist-dependent suppression of TNFSF14, NFKB2 and NFKBIZ was attenuated in NR1H3-KO THP-1 macrophages (Fig. 5K), supporting an NR1H3-dependent mechanism.

Cross-species integration highlights tree shrew as a phylogenetic bridge

Given the conserved role of NR1H3 in human and tree shrew macrophages, we investigated whether the identified regulatory mechanism represents a broader feature preserved across vertebrate evolution. To minimise confounding effects of tissue microenvironment-driven specialization, we focused on splenic macrophages. We integrated publicly available spleen scRNA-seq data across seven vertebrate species (GEO accession GSE186158), including primates (human and macaque), a close relative of primates (tree shrew), and commonly used experimental models spanning diverse lineages (pigs, rats, mice, and zebrafish). After quality control, the integrated map consisted of 88,106 cells (approximately 10,000 cells for each species), and was annotated into four major immune lineages based on canonical markers (Fig. 6A-B, Supplementary Fig. 18A). To compare cross‑species macrophage programs, we generated pseudobulk profiles by aggregating splenic macrophages within each species to reduce technical noise, and performed multidimensional scaling (MDS) on these profiles. The MDS projection broadly reflected phylogenetic relationships and positioned tree shrew between rodents and primates (Fig. 6C).

We first identified a core set of conserved macrophage genes broadly expressed across species, including CSF1R, NR1H3, CFP, CTSB, PSAP, LGMN, and ASAH1. CSF1R showed consistent expression across all species (Fig. 6D and E). NR1H3 expression varied across datasets, with low expression in mouse macrophages compared with higher expression in human and tree shrew macrophages (Fig. 6E). We next examined conservation of functional macrophage programs. Cytokine-related gene sets (IL1B, IL4R, IL10RA, TGFB1, IFNAR1, IFNGR1, TNFAIP8, TNFRSF14, and TNFRSF1B) were broadly conserved across species, with zebrafish showing a reduced set (TNFRSF14 and TGFB1). Chemokine gene (CCL4, CCL5, CXCR4, CXCL10, and CXCL16) showed closer similarity between tree shrew and human than between macaque and human in this dataset, differing mainly at CCL5. Complement-related genes (C1QA, C1QB, C1QBP, and C1QC) were largely conserved (Fig. 6E). To further explore the species-specific transcriptional features, we applied high-dimensional weighted gene coexpression network analysis (hdWGCNA)47 and identified 23 modules (Supplementary Fig. 18 B, C, Supplementary Data 6). We computed Pearson correlations between module genes and species labels, and designated the module with the highest correlation coefficient as the signature module for each species (Fig. 6F). Species-specific module genes were subjected to GO enrichment analysis using topGO (v2.38.1) (Fig. 6G and Supplementary Fig. 18 D). Although each species exhibited distinct functional characteristics associated with its macrophage-specific modules, many enriched terms converged on immune development and inflammatory pathways. Together, these results highlight both conservation and divergence of splenic macrophage programmes across vertebrates and support the tree shrew as a phylogenetically intermediate model for studying human-relevant immunological mechanisms.

Methods
Ethics statement

For the multi-organ single-cell atlas, tissues were collected from three adult tree shrews (two males and one female; 8 months old; 120.5 ± 4.3 g). For the EBV infection experiments, twelve tree shrews (8 months old; 123.1 ± 7.6 g) were randomly assigned to the 3 dpi and 7 dpi groups (n = 6 per group; three males and three females per group), and tissues were collected at the indicated time points. All tree shrews were housed separately in a semi-barrier facility at 23°C–25°C and 40%–50% humidity under a 12 h light/12 h dark cycle, with ad libitum access to standard chow and water60. All animal procedures were approved by the Animal Ethics Review Committee of Guangxi Medical University (approval no. 202101025) and were conducted in accordance with the Guiding Principles for the Use and Care of Experimental Animals issued by the Ministry of Science and Technology of China. Animals were euthanized by overdose of sodium pentobarbital followed by cervical dislocation prior to tissue collection.

Human peripheral blood samples were obtained from three healthy adult volunteers (two males and one female; 23-25 years old) for PBMC isolation. All participants provided written informed consent, and all uses of human material were approved by the Ethics Review Committee of the First Affiliated Hospital of Guangxi Medical University (approval no. 2025-E0760). The study was conducted in accordance with the Declaration of Helsinki and relevant local regulations.

Tissue dissociation and sample preparation for single-nucleus/cell library generation

Tissues were collected from 12 different samples: the spleen, thymus, tonsil, lymph nodes, nasopharynx, marrow, PBMCs, lung, liver, heart, kidney and stomach. Parenchymal tissues (except liver and heart) were placed in prechilled Hanks’ balanced salt solution (HBSS, Thermo Fisher Scientific) immediately after dissection. Samples were sliced into 3-4 mm3 pieces and digested into a cell suspension using digestion buffer containing 0.1 mg/mL DNase I (Sigma‒Aldrich) and 1 mg/mL collagenase II (Thermo Fisher Scientific). Cells were passed through a 70 μm strainer to remove the undigested tissue and centrifuged at 300 ×g for 5 min at 4 °C. Pelleted cells were resuspended in red blood cell lysis buffer (Solarbio) and incubated for 1 min at 4 °C. After washing twice with DPBS (Gibco), cell pellets were resuspended in sorting buffer (DPBS supplemented with 1% FBS). To mitigate dissociation-associated bias for large or fragile cells61, snRNA-seq was performed on liver and heart. Single-cell suspensions were captured in 10 × Genomics Chromium Single-Cell 5’ Solution, and single nuclei were captured in 10 × Genomics Chromium Single-cell 3’ Solution. The scRNA-seq/snRNA-seq library was generated by Beijing Novogene Company and prepared following the manufacturer’s protocol.

Raw data processing

Raw reads were aligned to the tree shrew reference genome (TupChi_1.0; https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000334495.1/) using the Cell Ranger v6.1.0 pipeline (10 × Genomics) with default parameters, followed by barcode assignment and unique molecular identifier (UMI) counting. Filtered count matrices were imported into Seurat (v4.0.0)62 for downstream processing. Potential doublets were detected and filtered using DoubletFinder63 with default parameters. The cells/nuclei were filtered on the basis of a minimal expression threshold of 500 genes and 1000 UMIs, as well as genes that were expressed in at least five cells/nuclei. To remove outliers, we discarded cells in the bottom 1% and top 1% of the gene and UMI distributions per library, and removed cells/nuclei with more than 30% mitochondrial counts. Filtered data were normalized with Seurat’s NormalizeData function and highly variable features (n = 2000) were selected using the “vst” method. Datasets across organs and modalities (scRNA‑seq/snRNA‑seq) were integrated using Seurat’s RPCA anchor‑based workflow. We computed 20 principal components (“RunPCA”) to mitigate batch effects, and generated a UMAP embedding for visualization. Clustering was performed using shared nearest neighbour graph clustering (Louvain community detection-based method) with “FindClusters” at resolution of 1.2. Cell types were annotated in two stages. First, major cell classes (immune, endothelial, epithelial, fibroblast, and parenchymal compartments) were assigned using canonical markers and cluster-defining genes identified with Seurat’s FindAllMarkers (default parameters). Second, immune‑cell subtypes were assigned using CellTypist (model: “Immune_All_High.pkl”; default settings) after mapping tree shrew genes to their human orthologs. Final labels were derived by a cluster‑level majority‑voting rule subject to a pre‑specified label‑probability threshold; cells not meeting this threshold were retained as “unassigned” and subsequently curated using canonical markers. Non‑immune subtypes were annotated using curated tissue‑specific marker panels from the literature and cross‑validated against the “FindAllMarkers”. Representative marker profiles for the 39 annotated cell types are shown in Supplementary Fig. 3.

Gene regulatory network construction

A list of tree shrew transcription factors was downloaded from AnimalTFDB3.064. We applied SCENIC (v0.10.0) to infer regulatory networks. TF-target pairs with weight values greater than 0.01 were retained to define coexpression modules. AUCell scores were computed to quantify regulon activity in each cell. The final cell type- and tissue-specific gene regulatory networks were visualized in Cytoscape 3.7.2 (https://github.com/cytoscape/cytoscape).

Cell–cell interaction analysis

Interactions among different cell types were predicted with the CellChat package (v1.1.3) based on curated ligand‒receptor pairs. Tree shrew genes were mapped to human orthologues (reciprocal best-hit BLASTP; Supplementary Data 7) to enable use of the human database. CellChatDB.human was used as the interaction database. The communication probabilities were calculated with the “computeCommunProb”, enriched ligand–receptor pairs were extracted using “extractEnrichedLR”. In addition, global communication patterns and signal networks were analysed using the built-in nonnegative matrix factorisation framework. Interactions with P value < 0.05 were considered significance.

Cross‑species comparisons

Public single-cell datasets for human, macaque, pigs, rats, mice, and zebrafish were obtained from the GEO database (accession number GSE186158) and processed using the same QC and normalization strategy prior to integration. For cross‑species analysis we restricted integration to spleen to minimize tissue‑driven heterogeneity and maintain consistency with the rest of the study. Orthologs were defined by reciprocal best‑hit BLAST of translated peptides against the human proteome; only one‑to‑one orthologs shared by all species were retained and gene identifiers were harmonized to human symbols (Supplementary Data 7). For each species, count matrices were processed in Seurat with SCTransform normalization, then performed integration using Seurat’s canonical correlation analysis (CCA) with default parameters. For between‑species transcriptomic similarity, we computed per‑species macrophage pseudobulk profiles, calculated Euclidean distances, and applied classical multidimensional scaling (MDS) to obtain the first two dimensions shown in the MDS plot46.

Statistical analysis

Data are presented as mean ± SD unless otherwise stated. Two-sided statistical tests were used throughout unless indicated. For two group comparisons, unpaired two-tailed Student’s t-test was used for approximately normally distributed data; otherwise the Wilcoxon rank-sum test was used. For comparisons involving more than two groups, one-way ANOVA with Tukey’s multiple-comparisons test (or Kruskal–Wallis with Dunn’s test for non-parametric data) was applied. Where multiple hypotheses were tested (e.g., differential expression, enrichment analyses), P values were adjusted using the Benjamini–Hochberg procedure, and FDR < 0.05 was considered significant unless stated otherwise.

Data availability

The scRNA-seq, bulk RNA-seq and Cut&Tag data generated in this study have been deposited in the NCBI Sequence Read Archive under BioProject accession PRJNA1051046 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1051046). Human splenic scRNA-seq was obtained from the Tissue Immune Cell Atlas portal (https://www.tissueimmunecellatlas.org/#datasets). Public datasets used for cross-species analyses were obtained from GEO (accession GSE186158). Source data are provided with this paper. Primer sequences and antibodies used in this study are provided in Supplementary Table 1 and Supplementary Table 2, respectively.

Code availability

The custom code used in this study is available at https://github.com/xialabtools/Tree_shrew_scRNA_code. This repository has been archived at Zenodo (https://doi.org/10.5281/zenodo.17609696).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript the authors analysed by single cell high-throughput RNA-seq several tissues from a small rodents (T.belangeri). The animal model was chosen because of a relatively higher similarity to humans with respect to mice and rats. Despite the inherently interesting findings regarding the tree shrew's immune system and its potential as a more representative model for human immunity, we have significant concerns regarding the reproducibility and the comprehensive functional validation of the study. The identification of human-like immune features, such as the tonsillar structure and the conserved NR1H3+ macrophage subset, is compelling, and the proposed role of these macrophages in EBV infection is particularly noteworthy.\n\nHowever, few critical issues must be addressed to solidify the study’s impact.\n\nGeneral comments:\nInsufficient Bioinformatics Pipeline Details and Reproducibility\nOur primary concern lies with the insufficient description and accessibility of the bioinformatics pipeline. While the provided \"Raw data processing\" section offers a good start by detailing some steps, it lacks the comprehensive transparency needed for true reproducibility, especially for complex single-cell analyses. The paper outlines exciting discoveries, but without a complete understanding and direct access to the computational methods used, it's impossible for other researchers to independently verify or build upon these findings.\n\nSpecifically, while the authors have listed some software and parameters, they need to provide a far more thorough and transparent account of their bioinformatics approach. This includes, but is not limited to:\nCustom Scripting and Workflows: The M&M section mentions specific packages and functions, but often single-cell analysis involves custom scripts to integrate these tools, perform specific filtering steps, or generate custom visualizations. These scripts, and the order in which they are executed, are crucial for replication.\n\nCell Type Annotation Strategy: A total of 39 cell types were identified in this paper, but the M&M section doesn't detail how these cell types were assigned. Was it based on known marker gene expression? Were external reference atlases used for annotation? This is a critical step that needs clear explanation.\n\nSpecifics of Cross-Species Comparison for CellChat: While it states \"homologous gene mapping from T. belangeri to the human genome was first performed to identify the corresponding genes in humans,\" the method or tool used for this homologous mapping (e.g., BLAST, Ensembl Compara, custom ortholog tables) is missing. This is vital for the validity of the cross-species comparison in CellChat.\n\nVersion Control and Environment Management: Simply listing software versions is helpful, but the exact computational environment (e.g., operating system, specific R/Python versions, library dependencies) in which the analysis was performed is essential for ensuring that code runs without dependency conflicts or unexpected behavior.\n\nCritically, to fully address the issue of reproducibility, I strongly recommend that the authors provide a link to a publicly accessible GitHub repository containing the complete analysis code. This repository should ideally be structured to allow for step-by-step execution of the analysis, from raw data processing to the generation of figures. Even better, the provision of a Docker container encapsulating the entire computational environment and analysis pipeline would be the gold standard.\nThis would allow any researcher to precisely replicate the computational steps and ensure the robustness of the reported findings, without encountering issues related to software installations or versioning.\n\n2. Lack of Functional Validation for NR1H3 Activation\nWhile the provided results presented in the paragraph “Conservation of NR1H3-mediated immunosuppression in human macrophages” demonstrate a strong correlation between NR1H3 activation and the downregulation of pro-inflammatory and non-canonical NF-KB pathway genes in THP-1 cells, it primarily focuses on transcriptional changes and direct binding. This is a valuable step, but it doesn't definitively prove the functional consequence of NR1H3 activation in mediating immunosuppression at a cellular or immunological level.\n\nTo demonstrate the \"true functionality\" of NR1H3 activation in immunosuppression, especially in the context of human macrophages, the following experiments would be helpful\na. Functional Assays on Macrophage Pro-inflammatory Responses:\nCytokine and chemokine Production Assays:\nMeasure secreted pro-inflammatory cytokines (using ELISA or multiplex cytokine arrays) upon stimulation with PAMPs/DAMPs from challenge NR1H3-OE, NR1H3-KO, and WT THP-1 cells\nNR1H3 activation should lead to a reduction in these cytokines compared to WT or NR1H3-KO cells.\nPhagocytosis Assay: While not directly linked to inflammation, some immunosuppressive macrophages might have altered phagocytic capabilities. This could be an exploratory experiment.\n\nb. Immune cell crosstalk assays:\nCo-culture Experiments with T cells (or other immune cells) can also give useful insight and functional readout. Therefore, a T cell Proliferation Assay can be performed to depict the biological consequences of NR1H3 activation. Likewise, the authors can investigate the T cell Cytokine Production and cell Differentiation\n\nc. In Vivo Relevance/Translational Experiments (if feasible):\nEx vivo stimulation of human primary macrophages: If possible, obtain human primary macrophages (e.g., from PBMCs) and activate NR1H3 pharmacologically (with T0901317 or GW3965). Then, repeat the cytokine/chemokine production assays after inflammatory stimulation. This would bridge the gap from THP-1 cells to primary human\nBy incorporating these functional assays, the authors would move beyond solely demonstrating transcriptional changes and direct gene binding to definitively prove the \"true functionality\" of NR1H3 activation in mediating immunosuppression in human macrophages. This would significantly strengthen their claim about the conserved anti-inflammatory immunoregulatory function of NR1H3.\n\nSpecific comments:\n1. The observed-to-expected cell ratios (Ro/e) presented in Fig 1 and Fig3D is not clear. What is the expected cell number? From which dataset does this derive? Wasn’t it worth to compare different tissues across themselves? As it is in this current version, the plot is not meaningful\n2. Fig 5A/B. The human/T.belangeri macrophage similarity is unclear and not convincing. What exactly is represented in Fig 5B? What is each dot? It is indeed difficult to state the “similarity”similarity”” between two cell types in different organs, unless many cells are compared and the two in exams are more close to each other with respect to the distance between other cells.\n\nMinor comments:\nIn the text, the experimental strategy is not clearly presented. The reader expects to have 36 samples, 12 tissues from three animals, but only 24 samples are actually there (3 for immune organs, one for other and one for pharyngeal tonsil). This must be better addressed at the beginning of the results.\nThe differences between gender is not analyzed in any of the main or supplementary figures. As there might be gender differences it is worth to analyze it and, especially for tissues where one single organ was processed, it is worth mentioning from which animal they were processed.\nFig 3e, the GO dot plot should be on the right of the GO terms, otherwise it seems the dots refers to the genes on the left\nLine 102 the goal is not to alleviate the financial constraints, so “and” should be replaced with “while”\nLine 158 fig 2e should be 2d\nLine 159 2f should be 2e\nFig 2C unit metrics is missing from scale bar, the same applies to most figures (2D, 3A, 3E, 4B, 4F, 4H and so on)\nThe volcano plot legend of Fig5E is not clear. It is more common to highlight the significantly upregulated genes (top right) and the significantly downregulated ones (top left)\nSupplementary figures are not mentioned in order in the text (i.e fig S11 comes after Fig S8)\nIn Fig6G human is misspelled\nFig 6C, this figure is not enough “per se” to state that T. belangeri is more close to humans than other species, as this plot might be dominated by strong species specific technical batch effect, rather than biological effect. Moreover, mice and rats seem to be closer to pigs than to themselves from this plot, which is counterintuitive and highlights possible misinterpretation or inadequate analysis.", "reviewer_2": "This study by Tang and colleagues presents a comprehensive single-cell immune cell atlas of Tupaia belangeri and provides compelling evidence that the tree shrew is a promising model organism for immunological research. By integrating multi-tissue profiling, EBV infection modeling, and cross-species comparisons, the authors effectively demonstrate the translational potential of Tupaia, particularly in contexts where conventional rodent models are limited.\n\nImportance:\nThe study is meaningful in that it establishes a foundational resource for using Tupaia belangeri as a non-primate animal model. Its evolutionary proximity to humans and permissiveness to human-specific pathogens make it particularly valuable for translational studies, especially in the field of viral immunology.\n\nSpecific comments:\n1. Data processing rigor:\nSeveral aspects of data processing and integration require clearer explanation. For instance, how doublets were removed, what specific QC metrics were applied, and how batch effects were corrected are not sufficiently described. Given the inclusion of multiple tissues and time points, this methodological clarity is essential.\nMore specifically, the localization of T cells in Supplementary Fig. 1 appears questionable, suggesting that doublet removal may not have been adequately performed. Epithelial cell annotations also appear inconsistent across figures. I recommend the authors include UMAP plots annotated by major cell types (e.g., T cells, B cells, epithelial cells, fibroblasts) to improve transparency.\nFurthermore, batch correction and integration should be benchmarked using multiple methods (e.g., Harmony, scVI) rather than relying on a single approach.\n\n2. Depth of analysis in the EBV model:\nThe EBV infection model is an important part of the study, but the analysis feels relatively shallow. For example, was EBV transcript detected in B cells or other cell types? Basic analyses such as pseudotime trajectory, cell type-specific differential gene expression (DEG), and pathway enrichment should be provided before diving into cell–cell communication analysis.\nAdditionally, recent tools like Milo could be combined with improved integration methods to track graph-level changes in cell proportions.\n\nMinor comment: In Supplementary Fig. 13A, the temporal order should be displayed as Control → Day 3 → Day 7. Also, the terms “Pos-macro” and “Neg-macro” are not clearly defined and are inconsistent with previous labeling.\n\n3. Cross-species comparison methodology:\nThe cross-species integration is one of the more exciting aspects of this paper, but its current implementation is too preliminary. Key technical steps—such as ortholog selection, normalization across species, and gene filtering—should be explained in more detail. Methods like SATURN or scGLUE could be applied here to improve both the interpretability and robustness of the analysis.\nIn particular, the MDS plot is difficult to interpret, and it is unclear what the axes represent. If this plot is intended to convey evolutionary distance or transcriptional similarity, it should be supplemented with more formal statistical support.\nAdditionally, the analysis of HLA loci, which are known to be highly divergent between species, needs to be handled with extra caution. Ortholog mapping for HLA genes should be genome-informed and well-documented.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "This new study describes a single cell immune atlas for the tree shrew (Tupaia belangeri). Tree shrews are an emerging small animal model that are evolutionarily closer to primates than mice. As such, tree shrews represent a useful, viable alternative model for infectious disease and aging studies that may be otherwise challenging to carryout in murine or non-human primate model systems. The investigators profiled single cells from 12 different tissues of 3 adult tree shrews, identifying 39 different cell types. Notably, bulk RNAseq studies have been performed in spleen, PBMCs, and other tissues from the tree shrew, and single cell transcriptome atlases have been recently published for brain and neurological tissues. A strength of this study is that it represents the first tree shrew immune-focused single cell atlas which can provide an important platform for understanding systemic transcriptional responses to infectious diseases in future studies. Aside from the substantial efforts that went into generating this resource, there are some key weaknesses and major concerns with the overall conclusions derived from the interpretation of the EBV infection data.\n\nMajor Concerns:\n1. Atlas reveals that T.belangeri possess well-developed tonsillar structures analogous to human tonsils. These tonsillar anatomical structures within the oropharynx were previously reported for the tree shrew—and thus, it is expected that these tissues will harbor lymphoid cell types commonly observed in primate tonsils. However, in figure 1, mature and activated B cell subsets are not detected in the pharyngeal tonsils, which is unusual. Could this be due to sampling issues or issues with adequate sequencing coverage? As the tonsils are a significant focus of the paper (Fig 2), more work needs to be put into describing these cells. Fig. 2B needs to be described in detail—UMAP represents cells in which tonsils? There were only ~2700 cells sequenced for pharyngeal—were palantine and pharyngeal grouped together in this figure? Are the cell subsets and relative distributions of representative cell types in these tissues comparable to what is reported in human and nonhuman primate samples?\n\n2. A major caveat is that the EBV infection model does not mimic natural oral uptake and natural viral transmission which is observed in other NHP models with LCVs (i.e. rhesus and cynomolgus macaque models). Here, the route of EBV administration is intravenous into the femoral vein, thus any relevance to palantine/pharyngeal tonsil responses—which the authors argue is an important component of the tree shrew model (Fig. 2)—is unclear as this requires trafficking of the virus to the oral/throat cavity (in previous publications with this tree shrew model, EBV shedding occurs at 7 dpi). It is unclear why tonsil cells were not sequenced from the infected animals as this would highlight the importance and biological relevance of these structures for EBV infection. The I.V. infection route and also the high titer virus administered is akin to what is utilized in mouse models, which does not recapitulate natural EBV infection in humans and this needs to be acknowledged. As such, the immune responses (i.e. bypassing mucosal immunity in the oral cavity and gut that occurs with oral challenge) are expected to be skewed.\n\n3. The sample size for EBV experiments in Fig 4 is extremely small. Spleen cells are sequenced from only two animals that were infected with EBV and harvested at different time points, limiting the statistical power and generalizability of findings that NR1H3(LXRa+) macrophages are indeed altered by EBV infection. It is basically a pilot experiment. In an attempt to confirm scRNAseq findings, immunofluorescence and qRT-PCR experiments were performed on samples from 4 animals in Fig 4. No differences in NR1H3+ or IRF8+ macrophage numbers were detected by IFA, although the NR1H3 RNA level does decrease, which hardly confirms the scRNAseq findings. Given that data describing NR1H3 cell reprogramming is from n=1, additional evidence should be provided to support the major conclusions that NRH13+ macrophages dominate the inflammatory responses observed in EBV infection. For example, does depletion of these cells impact acute EBV infection? Are there additional signs of inflammatory markers that track with these macrophage subsets?\n\n4. Depletion of liver X receptors (NR1H3 gene encodes LXRa) in mice leads to enhanced replication of a murine g-herpesvirus (MHV68) that is often used as a model for EBV infection (PMC6050960). Such prior studies indicate NRH13-mediated mechanisms are indeed conserved in mice.\n\n5. The study largely focuses on describing B cells and macrophage subset, yet T cell responses are a major contributor to control of infectious agents that have been tested in tree shrew models (i.e. Zika, hepatitis virus). T cell populations should be described in detail.\n\nOther minor comments:\n1. Introduction, line 52: smaller NHPs such as marmosets and tamarins have been historically used for EBV infection studies and should be mentioned in this section.\n2. Line 74, first single-cell atlas of T. belangeri and line 435: This is incorrectly stated as there are at least three other publications describing brain/neurological tissue-focused single cell atlases for the tree shrew."}, "author_response": "We appreciate the reviewers for their insightful comments and suggestions, which helped to improve the quality of our work. We have made substantial revision to address all of your comments. Point-by-point responses to the reviewers’ comments and concerns are presented below.\n\nReviewers’ Comments to Author\nReviewer #1:\nRemarks to the Authors\nIn this manuscript the authors analysed by single cell high-throughput RNA-seq several tissues from a small rodents (T.belangeri). The animal model was chosen because of a relatively higher similarity to humans with respect to mice and rats. Despite the inherently interesting findings regarding the tree shrew's immune system and its potential as a more representative model for human immunity, we have significant concerns regarding the reproducibility and the comprehensive functional validation of the study. The identification of human-like immune features, such as the tonsillar structure and the conserved NR1H3+ macrophage subset, is compelling, and the proposed role of these macrophages in EBV infection is particularly noteworthy.\nHowever, few critical issues must be addressed to solidify the study’s impact.\nResponse: Thanks for your comments concerning our paper. Those comments are all valuable and very helpful for revising and improving our paper, as well as the important guiding significance to our research. We have studied comments carefully and have made modification which we hope meet with approval. Revised portion are marked in the paper.\nGeneral comments:\n1. Insufficient Bioinformatics Pipeline Details and Reproducibility: Our primary concern lies with the insufficient description and accessibility of the bioinformatics pipeline. While the provided \"Raw data processing\" section offers a good start by detailing some steps, it lacks the comprehensive transparency needed for true reproducibility, especially for complex single-cell analyses. The paper outlines exciting discoveries, but without a complete understanding and direct access to the computational methods used, it's impossible for other researchers to independently verify or build upon these findings. Specifically, while the authors have listed some software and parameters, they need to provide a far more thorough and transparent account of their bioinformatics approach. This includes, but is not limited to: (a) Custom Scripting and Workflows: The M&M section mentions specific packages and functions, but often single-cell analysis involves custom scripts to integrate these tools, perform specific filtering steps, or generate custom visualizations. These scripts, and the order in which they are executed, are crucial for replication.\nResponse: We thank the reviewer for emphasizing the importance of transparency and reproducibility. We fully agree and have now released a public GitHub repository (https://github.com/xialabtools/Tree_shrew_scRNA_code) together with a long‑term archive on Zenodo (https://doi.org/10.5281/zenodo.17609696). The repository contains an end‑to‑end, executable workflow from raw data to all figures, including raw‑data processing and QC, doublet detection, data integration and annotation, cross‑species analyses, and per‑figure reproduction scripts. Each step specifies the exact software used and key parameters, and the order of execution is documented to enable independent replication.\nChanges in the manuscript: We have added a “Code availability” statement at Line 734: The custom code used in this study is available at https://github.com/xialabtools/Tree_shrew_scRNA_code. This repository has been archived at Zenodo (https://doi.org/10.5281/zenodo.17609696).”\n(b) Cell Type Annotation Strategy: A total of 39 cell types were identified in this paper, but the M&M section doesn't detail how these cell types were assigned. Was it based on known marker gene expression? Were external reference atlases used for annotation? This is a critical step that needs clear explanation.\nResponse: We thank the reviewer for highlighting the need to clarify our annotation strategy. We implemented a two-stage procedure: (i) We first separated major classes using canonical markers (e.g., immune vs. non-immune; endothelial, epithelial, fibroblast, parenchymal) at the cluster level, guided by cluster-defining genes identified with Seurat’s FindAllMarkers (default parameters). This step leverages the marker framework already described in the manuscript and summarized across the 39 cell types. (ii) Immune compartment: we applied CellTypist1 (model “Immune_All_High.pkl”, default settings) after mapping tree shrew genes to human orthologs, and adopted cluster‑level majority voting with a label‑probability threshold to retain only high‑confidence labels. Non‑immune compartment: we used curated organ‑specific canonical markers from the literature to assign labels and cross‑checked them with the FindAllMarkers outputs. Ambiguous clusters were flagged for manual curation. For transparency, the full per‑cluster marker tables, annotation scripts, and execution order are provided in our public repository/Zenodo snapshot (see Code Availability). The final set of 39 annotated cell types and representative marker patterns is summarized in Supplementary Fig. 3, which supports the validity of our labels.\nChanges in the manuscript: The following text has been added Raw data Line 522: “Cell types were annotated in two stages. First, major cell classes (immune, endothelial, epithelial, fibroblast, and parenchymal compartments) were assigned at the cluster level on the basis of canonical markers, informed by cluster‑defining genes identified with Seurat’s FindAllMarkers (default parameters). Second, immune‑cell subtypes were assigned using CellTypist (model: “Immune_All_High.pkl”; default settings) after mapping T. belangeri genes to their human orthologs. Final labels were derived by a cluster‑level majority‑voting rule subject to a pre‑specified label‑probability threshold. Cells not meeting this threshold were retained as “unassigned” and subsequently curated on the basis of canonical markers. Non‑immune subtypes were annotated using curated organ‑specific marker panels from the literature and cross‑validated against the FindAllMarkers outputs. Representative marker profiles for the 39 annotated cell types are shown in Supplementary Fig. 3.”\n(c) Specifics of Cross-Species Comparison for CellChat: While it states \"homologous gene mapping from T. belangeri to the human genome was first performed to identify the corresponding genes in humans,\" the method or tool used for this homologous mapping (e.g., BLAST, Ensembl Compara, custom ortholog tables) is missing. This is vital for the validity of the cross-species comparison in CellChat.\nResponse: We thank the reviewer for this valuable comment. For the CellChat analyses, because the ground truth database was curated from the human literature and insufficient data exist to curate ligand-receptor pairs for each species. T.belangeri transcriptomes were converted into human-equivalent transcriptomes by mapping of genes to their top human orthologue as detected by reciprocal best-hit BLASTP search. This transformed expression matrix was used only for cell communication analysis: dimensionality reduction, cluster identification, differential expression testing, and marker gene identification and plotting were conducted on unadulterated species-relevant gene sets. This workflow is consistent with the practice for cross‑species single‑cell comparisons and has been used in similar contexts (e.g., Tarashansky et al., 20212). A complete ortholog table is provided in the Supplementary Table7.\nChanges in the manuscript: In Materials and Methods (section “Cell-cell interaction analysis”, Line 537), the following text has been added”: “To perform CellChat analysis, T. belangeri transcriptomes were converted into human-equivalent transcriptomes by mapping of genes to their top human orthologue as detected by reciprocal best-hit BLASTP search (Supplementary Table7).”\n(d) Version Control and Environment Management: Simply listing software versions is helpful, but the exact computational environment (e.g., operating system, specific R/Python versions, library dependencies) in which the analysis was performed is essential for ensuring that code runs without dependency conflicts or unexpected behavior.\nResponse: We thank the reviewer for stressing the importance of a well‑specified computing environment. Beyond the detailed Methods in the manuscript, we provide a public GitHub repository that contains the complete analysis code from raw matrices to figure scripts, together with an R sessionInfo summary of the software stack to document the exact package versions and platform used. These resources should enable independent researchers to reproduce the analyses as reported.\n(e) Critically, to fully address the issue of reproducibility, I strongly recommend that the authors provide a link to a publicly accessible GitHub repository containing the complete analysis code. This repository should ideally be structured to allow for step-by-step execution of the analysis, from raw data processing to the generation of figures. Even better, the provision of a Docker container encapsulating the entire computational environment and analysis pipeline would be the gold standard. This would allow any researcher to precisely replicate the computational steps and ensure the robustness of the reported findings, without encountering issues related to software installations or versioning.\nResponse: We agree that transparent access to the full analysis is essential. The step‑wise repository (linked in the Code availability statement) enables execution from raw data processing to figure generation, and a versioned release has been archived with a permanent DOI to ensure long‑term accessibility. While container images are not required for our HPC setting, the combination of (i) explicit software and parameterization in the Methods, (ii) per‑step scripts in the repository, and (iii) an R sessionInfo summary provides a complete and auditable path for reproduction. The corresponding author will gladly assist with any additional technical details if needed.\nChanges in the manuscript: We have added at Line 734: “Code availability: The custom code used in this study is available at https://github.com/xialabtools/Tree_shrew_scRNA_code. This repository has been archived at Zenodo (https://doi.org/10.5281/zenodo.17609696).”\n2. Lack of Functional Validation for NR1H3 Activation\nWhile the provided results presented in the paragraph “Conservation of NR1H3-mediated immunosuppression in human macrophages” demonstrate a strong correlation between NR1H3 activation and the downregulation of pro-inflammatory and non-canonical NF-KB pathway genes in THP-1 cells, it primarily focuses on transcriptional changes and direct binding. This is a valuable step, but it doesn't definitively prove the functional consequence of NR1H3 activation in mediating immunosuppression at a cellular or immunological level. To demonstrate the “true functionality” of NR1H3 activation in immunosuppression, especially in the context of human macrophages, the following experiments would be helpful.\n(a) Functional Assays on Macrophage Pro-inflammatory Responses: (1) Cytokine and chemokine Production Assays: Measure secreted pro-inflammatory cytokines (using ELISA or multiplex cytokine arrays) upon stimulation with PAMPs/DAMPs from challenge NR1H3-OE, NR1H3-KO, and WT THP-1 cells. NR1H3 activation should lead to a reduction in these cytokines compared to WT or NR1H3-KO cells.\nResponse: Thank you for this constructive suggestion. We added a cytokine‑secretion assay to test whether NR1H3 activation suppresses macrophage inflammatory cytokines. Specifically, we challenged WT, NR1H3‑overexpressing (OE), and NR1H3‑knockout (KO) THP‑1–derived macrophages with a standardized PAMP stimulus and quantified secreted IL‑1β, IL‑6, and TNF‑α in culture supernatants by ELISA (Elabscience, China); absorbance was read at 450 nm on a BioTek microplate reader with on‑plate standards confirming good linearity (R² > 0.98). Across n=4 biological replicates per group and one‑way ANOVA with post‑hoc testing, NR1H3‑OE cells released significantly less IL‑1β, IL‑6, and TNF‑α than WT and KO (P < 0.01 for each cytokine). KO vs WT showed a significant increase for IL‑1β (P < 0.01) and no significant change for IL‑6/TNF‑α under this strong stimulus (likely a ceiling effect) (RFigure 1). These secretion data provide the requested functional readout and are consistent with our human‑cell mechanism showing NR1H3‑dependent repression of TNF and NF‑κB genes (NFKB2, NFKBIZ, MAP3K14) with direct NR1H3 binding, thereby linking molecular repression and reduced effector cytokines in human macrophages.\nRFigure 1: ELISA validation of NR1H3-mediated suppression of cytokine secretion in THP-1 macrophages. Left: ELISA plate maps (WT, NR1H3-KO, NR1H3-OE) with diluted standards. Middle: standard curves with linear regression and R² (IL-1β: 0.9925; IL-6: 0.9954; TNF-α: 0.9876). Right: cytokine quantification (mean ± SD, n=4); P < 0.01 indicated by **, “ns” denotes not significant (one-way ANOVA with post-hoc tests).\n(2) Phagocytosis Assay: While not directly linked to inflammation, some immunosuppressive macrophages might have altered phagocytic capabilities. This could be an exploratory experiment.\nResponse: To explore whether NR1H3 activation alters macrophage effector functions, we compared phagocytosis in WT vs NR1H3‑overexpressing (OE) THP‑1–derived macrophages. THP‑1 cells (5×105/mL) were differentiated with PMA 100 ng/mL for 24 h, then incubated 4 h with pHrodo™ Red E. coli bioparticles (Thermo Fisher; 1:100 v/v; reagent resuspended in PBS and sonicated 5 min). Phagocytosed particles (fluorescent in acidified phagosomes) were visualized by confocal microscopy and quantified by flow cytometry. NR1H3‑OE cells showed a lower pHrodo‑positive fraction than WT (representative gates: WT = 91.6%, OE = 64.4%), with a significant reduction in phagocytic cells (P < 0.01, unpaired two‑tailed t‑test; mean ± SD; independent experiments shown as dots)( RFigure 2).\nThese data suggest NR1H3 activation modestly dampens phagocytic uptake, complementing our molecular evidence that NR1H3 enforces an immunoregulatory program in human macrophages.\nRFigure 2: NR1H3 activation reduces phagocytosis of pHrodo‑labeled bacteria in THP‑1 macrophages. (A) Phase‑contrast/fluorescence overlays of WT and NR1H3‑OE cells exposed to pHrodo Red E. coli (red). (B) Confocal images (DAPI, blue; pHrodo, red). (C) Representative flow‑cytometry histograms showing gated pHrodo‑positive fractions (left: WT + E. coli; right: OE + E. coli) and summary quantification (right panel; mean ± SD). ** P < 0.01; unpaired two‑tailed t‑test. Scale bars: 100 µm.\n(b) Immune cell crosstalk assays:\nCo-culture Experiments with T cells (or other immune cells) can also give useful insight and functional readout. Therefore, a T cell Proliferation Assay can be performed to depict the biological consequences of NR1H3 activation. Likewise, the authors can investigate the T cell Cytokine Production and cell Differentiation\nResponse: To test whether NR1H3 activation in human macrophages functionally skews T‑cell responses, we isolated CD3⁺ T cells from healthy donors (written consent; Ethics 2025‑E0760) and stimulated them with PHA in the presence of conditioned supernatants from either WT‑THP‑1 or NR1H3‑overexpressing (OE) THP‑1 macrophages. IFN‑γ output was measured by ELISPOT; T‑cell effector polarization was assessed by intracellular IFN‑γ staining within CD8⁺ T cells using a standard Fc‑block, surface stain, fixation/permeabilization and intracellular staining on a CytoFLEX; and proliferation was quantified by CFSE dilution (RFigure 3A).\nResults: NR1H3‑OE macrophage supernatant reduced IFN‑γ spots compared with WT (P<0.05) (RFigure 3B), and decreased the fraction of IFN‑γ⁺ CD8⁺ T cells (54.8% to 28.2%, P<0.01) (RFigure 3C), while proliferation was unchanged (CFSE, not significant) (RFigure 3D). Thus, NR1H3 activation in macrophages dampens T‑cell effector cytokine production/polarization without nonspecific suppression of clonal expansion, providing a clear cellular‑level consequence that is consistent with our demonstrated NR1H3‑dependent repression of inflammation function in human macrophages.\nRFigure 3: NR1H3‑activated macrophages restrain T‑cell IFN‑γ responses without affecting proliferation. (A) Workflow schematic: donor PBMC isolation, CD3⁺ T‑cell negative selection, culture with WT or NR1H3‑OE THP‑1 conditioned supernatants during PHA stimulation, and ELISPOT/flow‑cytometry readouts. (B) IFN‑γ ELISPOT images (left) and quantification (right, mean ± SD; P<0.05). (C) Flow cytometry: representative plots showing IFN‑γ⁺ cells within CD8⁺ T cells (WT = 54.8%, OE = 28.2%); right, summary (mean ± SD); ** P<0.01. (D) CFSE‑based proliferation: overlaid histograms and quantification indicate no significant difference (n.s.).\n(c) In Vivo Relevance/Translational Experiments (if feasible):\nEx vivo stimulation of human primary macrophages: If possible, obtain human primary macrophages (e.g., from PBMCs) and activate NR1H3 pharmacologically (with T0901317 or GW3965). Then, repeat the cytokine/chemokine production assays after inflammatory stimulation. This would bridge the gap from THP-1 cells to primary human\nResponse: To extend our findings from THP‑1 cells to primary human cells, we isolated PBMCs from healthy donors with written consent and ethics approval 2025‑E0760. Total monocytes were purified by negative selection and differentiated into macrophages for 5 days with M‑CSF (50 ng/mL); adherent cells were then treated with the NR1H3 agonist T0901317 for 24 h, challenged with an inflammatory stimulus, and culture supernatants were assayed for IL‑1β, IL‑6, and TNF‑α by ELISA (Elabscience; OD 450 nm, BioTek) (RFigure 4A). Across independent donors, T0901317 significantly lowered all three cytokines compared with control macrophages (P < 0.01) (RFigure 4B). Together, these ex vivo human data demonstrate that pharmacologic NR1H3 activation suppresses primary macrophage inflammatory output, reinforcing the translational relevance of our THP‑1 and mechanistic results showing NR1H3‑dependent repression of TNF and non‑canonical NF‑κB programs in human macrophages.\nRFigure 4: NR1H3 activation suppresses cytokine output in primary macrophages. (A) Workflow: donor PBMC isolation; monocyte negative selection; macrophage differentiation with M‑CSF; NR1H3 activation with T0901317; inflammatory challenge; ELISA readouts and T‑cell assay. (B) ELISA of T0901317‑activated macrophages (mean ± SD; dots represent independent experiments; *P < 0.05, ** P < 0.01).\nBy incorporating these functional assays, the authors would move beyond solely demonstrating transcriptional changes and direct gene binding to definitively prove the “true functionality” of NR1H3 activation in mediating immunosuppression in human macrophages. This would significantly strengthen their claim about the conserved anti-inflammatory immunoregulatory function of NR1H3.\nResponse: We are grateful for this constructive recommendation, which substantially improved the manuscript. These additions go beyond transcriptional and binding evidence and provide direct cellular and ex vivo confirmation of the true functionality of NR1H3 activation in mediating immunosuppression in human macrophages, thereby strengthening our claim of a conserved anti‑inflammatory role for NR1H3. In response, we performed the requested functional assays and compiled the results as Supplementary Figure 16 (RFigure 5). We also revised the main text at Line 287 to read: “In functional assays, NR1H3 activation in THP‑1 macrophages reduced secretion of IL‑1β, IL‑6 and TNF‑α, modestly decreased uptake of bacterial particles, and diminished IFN‑γ production and effector polarization of primary CD8⁺ T cells while leaving T‑cell proliferation unchanged (Supplementary Fig. 18A–F). Extending these findings ex vivo, pharmacologic activation of NR1H3 with T0901317 in primary human monocyte‑derived macrophages similarly suppressed LPS‑induced IL‑1β, TNF‑α and IL‑6 release (Supplementary Fig. 18G,H), underscoring the translational relevance of NR1H3‑mediated immunoregulation.”\nRFigure 5: Functional of NR1H3 activation in human macrophages and translational ex vivo validation. (A) THP 1–derived macrophages (WT, NR1H3 KO and NR1H3 OE) were stimulated under a standardized pro inflammatory challenge; supernatants were collected and IL 1β, IL 6 and TNF α were quantified by ELISA (Elabscience; OD 450 nm, BioTek). NR1H3 OE released significantly less of all three cytokines than WT and KO, whereas KO showed an increase for IL 1β and no significant change for IL 6/TNF α under these conditions. (B) Phagocytosis assay using pHrodo™ Red E. coli bioparticles (Thermo Fisher). THP 1 cells were differentiated with PMA 100 ng/mL for 24 h, incubated with bioparticles (1:100, 4 h), and analyzed by flow cytometry; NR1H3 OE displayed a reduced pHrodo positive fraction compared with WT. (C) Schematic of the T cell crosstalk workflow: CD3⁺ T cells were negatively selected from healthy donor PBMCs (written consent; Ethics 2025 E0760), PHA stimulated, and cultured with conditioned supernatants from WT or NR1H3 OE THP 1 macrophages. (D) IFN γ ELISPOT showed fewer spots in T cells exposed to NR1H3 OE supernatant than to WT supernatant. (E) Flow cytometry with intracellular cytokine staining demonstrated a lower frequency of IFN γ⁺ cells within CD8⁺ T cells after exposure to NR1H3 OE supernatant. (F) CFSE dilution analysis indicated no difference in T cell proliferation between conditions. (G) Ex vivo validation in primary human monocyte derived macrophages: CD14⁺ monocytes were isolated by negative selection, differentiated for 5 days with M CSF (50 ng/mL) to adherent macrophages, treated with the NR1H3 agonist T0901317 for 24 h, then challenged with LPS as indicated. (H) ELISA of macrophage supernatants showed significant reductions in IL 1β, TNF α and IL 6 after T0901317 treatment compared with control. Data are mean ± SD; each dot represents an independent experiment or donor. *P<0.05; **P<0.01; ns, not significant.\nSpecific comments:\n1. The observed-to-expected cell ratios (Ro/e) presented in Fig 1 and Fig3D is not clear. What is the expected cell number? From which dataset does this derive? Wasn’t it worth to compare different tissues across themselves? As it is in this current version, the plot is not meaningful\nResponse: We thank the reviewer for raising this point. In our study, the observed‑to‑expected ratio (Ro/e) quantifies tissue preference for each cell type under a common global baseline. For cell type i in tissue t, the expected count is defined as E = pi × Nt, where pi is the global proportion of that cell type computed from our i,t integrated dataset of all tissues in this study, and Nt is the total number of profiled cells in tissue t; the observed count Oi,t is the number of cells of type i in t. The enrichment ratio is Ro/e = O / E ; values > 1 indicate enrichment and <1 indicate i,t i,t i,t depletion. Thus, Ro/e explicitly compares tissues against each other for each cell type using a common global baseline, which is what our heatmaps visualize (e.g., Fig. 1D and Fig. 3D; thresholds noted in the legends to highlight strong preference). This formulation follows prior single-cell work (e.g., Zhang et al., Nature. 2018)3 and is referenced in our manuscript; we now make the formula and data source explicit in the text and figure legends for clarity.\nChanges in the manuscript: We have revised the Methods (Statistical analysis) Line 554 to explain that expected counts are based on the global cell-type distribution: Expected count for cell type i in tissue t = (global proportion of i across all tissues) × (total cells in tissue t)41. Ro/e > 1 indicates enrichment, Ro/e <1 indicates depletion.\n2. Fig 5A/B. The human/T.belangeri macrophage similarity is unclear and not convincing. What exactly is represented in Fig 5B? What is each dot? It is indeed difficult to state the “similarity” between two cell types in different organs, unless many cells are compared and the two in exams are more close to each other with respect to the distance between other cells.\nResponse: We appreciate the opportunity to clarify this figure. Fig. 5A shows the NR1H3⁺ macrophage subset in the human spleen dataset (Domínguez‑Conde et al., 2022; Ref. 43), not cells from other organs. Fig. 5B quantitatively assesses cross‑species transcriptional similarity between splenic NR1H3⁺ macrophages in human and T. belangeri. Specifically, we (i) restrict both datasets to one‑to‑one orthologs, (ii) compute pseudobulk expression profiles for NR1H3⁺ macrophages in each species (averaging across all cells to obtain robust profiles), and (iii) calculate a Pearson correlation between the two pseudobulk profiles (r = 0.77, p < 0.01). Each dot in Fig. 5B corresponds to one ortholog‑matched gene (mean log‑normalized expression in human on one axis and in T. belangeri on the other). We have revised the legend of Fig. 5B accordingly to avoid any ambiguity. This procedure and the resulting correlation are also described in the Results.\nChanges in the manuscript: We have revised the Fig. 5 B legend Line 663: “Pearson correlation (r = 0.77, P < 0.01) between pseudobulk expression profiles of splenic NR1H3⁺ macrophages from human and T. belangeri after one‑to‑one ortholog mapping; each dot represents the mean log‑normalized expression of an orthologous gene (averaged across cells), and the axes show the mean log‑normalized expression in each species.”\nMinor comments:\n(a) In the text, the experimental strategy is not clearly presented. The reader expects to have 36 samples, 12 tissues from three animals, but only 24 samples are actually there (3 for immune organs, one for other and one for pharyngeal tonsil). This must be better addressed at the beginning of the results.\nResponse: We thank the reviewer for highlighting this ambiguity in the study design. We did not generate 36 independent libraries (12 tissues × 3 animals). Instead, we generated 24 single‑cell/single‑nucleus libraries across 12 organs from three adult animals, with heart and liver profiled by snRNA‑seq to preserve fragile/large cells. To maximize immunological insight under realistic resources, immune‑related organs were sampled across animals, whereas non‑immune parenchymal organs and the pharyngeal tonsil were profiled once as representative tissues. This design choice and the rationale are now stated explicitly at the start of the Results, with per‑tissue cell numbers and modalities summarized in Fig. 1B and Supplementary Table 1. We also direct readers to Supplementary Fig. 2, which shows concordant cell‑type distributions among the same/related immune tissues, supporting reproducibility despite uneven replication.\nChanges in the manuscript: At the beginning of the Results (Line 93), we added: “We generated 24 sc/snRNA‑seq libraries across 12 organs from three adult tree shrews; immune‑related organs were sampled across animals, whereas non‑immune parenchymal organs and the pharyngeal tonsil were profiled once due to study scope and tissue availability (see Fig. 1B and Supplementary Table 1).” We also added a brief note in Methods (“Tissue dissociation and sample preparation” and “Raw data processing”) clarifying scRNA‑seq vs snRNA‑seq usage (heart/liver) and why independent replicates were limited outside the immune system.\n(b) The differences between gender is not analyzed in any of the main or supplementary figures. As there might be gender differences it is worth to analyze it and, especially for tissues where one single organ was processed, it is worth mentioning from which animal they were processed.\nResponse: We appreciate the suggestion to clarify sex information and potential sex effects. Our study comprised three adult animals (two males and one female), as now stated in the Methods and summarized per tissue in Supplementary Table 1. Multiple immune tissues (e.g., pharyngeal tonsil, spleen, lymph node) were sampled from all three animals, so analyses in these organs include data from both sexes. Some non‑immune tissues (e.g., liver, heart) were profiled in a single donor; the donor’s sex is explicitly indicated in Supplementary Table 1 and the relevant figure legends. Given this imbalanced sample size (1 female vs 2 males), we are underpowered for formal sex‑stratified differential expression. We therefore did not pursue sex‑specific DE as a main analysis but now acknowledge this as a limitation. For transparency, we provide a new supplementary UMAP split by sex showing no sex‑driven segregation of the global embedding, which reduces concern that our principal conclusions are confounded by sex in this dataset (RFigure 6).\nChanges in the manuscript: We have added (Discussion Limitations Line 439): “We mention that sex-based differential expression was not explored due to the small number of animals of each sex, and we acknowledge this as a limitation of the current dataset.”\nRFigure 6: Sex-stratified UMAP of the integrated tree shrew atlas.\n(c) Fig 3e, the GO dot plot should be on the right of the GO terms, otherwise it seems the dots refers to the genes on the left\nResponse: We have revised Fig. 3E to remove any ambiguity. The marker‑gene block has been moved to the far left of panel E and is visually separated from the enrichment results; the GO bubble plot is now aligned directly with the GO‑term list on the right, and the legend clarifies that bubble size = gene count, bubble color = gene ratio, and x‑axis = −log10(adjusted P). This makes it clear that bubbles correspond to GO terms; the underlying analysis is unchanged. If the reviewer prefers the alternative ordering (placing the bubbles to the right of the GO‑term labels), we can provide that layout as well.\n(d) Line 102 the goal is not to alleviate the financial constraints, so “and” should be replaced with “while”\nResponse: Thank you. We have replaced “and” with “while” as suggested.\n(e) Line 158 fig 2e should be 2d\nResponse: Thank you. We have corrected the “fig 2e” with “fig 2d”.\n(f) Line 159 2f should be 2e\nResponse: Thank you. We have corrected “fig 2f” with “fig 2e”.\n(g) Fig 2C unit metrics is missing from scale bar, the same applies to most figures (2D, 3A, 3E, 4B, 4F, 4H and so on)\nResponse: Thank you for pointing this out. We audited all cited panels and added missing labels/units to make the scales explicit. Specifically:\n1) Figs. 2C, 4B, 4H (dot plots): colour bar now labelled “log‑normalized expression”.\n2) Fig. 3A (feature/UMAP): colour bar labelled “log‑normalized expression.”\n3) Fig. 2D (KEGG heatmap): colour bar labelled “−log10(adjusted P value).”\n4) Fig. 3E: heatmap labelled “Row z‑score,” GO bubble panel explicitly indicates gene ratio (colour), gene count (size), and x‑axis = −log10(adjusted P).\nChanges in the manuscript: We have added the labels/units above to the corresponding panels and updated the figure legends to reflect them, including an explicit note in Fig. 4F that pseudotime has no units.\n(h) The volcano plot legend of Fig5E is not clear. It is more common to highlight the significantly upregulated genes (top right) and the significantly downregulated ones (top left)\nResponse: We agree that the legend could be clearer. We have revised Fig. 5E to follow the common convention: x‑axis = log2 fold change (NR1H3‑OE vs WT; positive = up in NR1H3‑OE), y‑axis = −log10(P). Points are colored red (significantly upregulated; P < 0.05 and log2FC ≥ 1), blue (significantly downregulated; P < 0.05 and log2FC ≤ -1), and grey (not significant). Vertical dashed lines mark |log2FC| = 1 and the horizontal dashed line marks P = 0.05. Representative genes in both directions remain annotated. The underlying statistics are unchanged.\nChanges in the manuscript: Updated Fig. 5E and legend accordingly.\n(i) Supplementary figures are not mentioned in order in the text (i.e fig S11 comes after Fig S8)\nResponse: We appreciate the careful reading. We have corrected the citation order. At Line 134, immediately after citing Supplementary Figs. S6–S7, we now cite Supplementary Figs. S8–S9. In addition, the previous miscitation of “Supplementary Fig. S8” has been corrected to “Supplementary Fig. S10”.\n(j) In Fig6G human is misspelled\nResponse: We have corrected the spelling of “human” in Fig. 6G.\n(k) Fig 6C, this figure is not enough “per se” to state that T. belangeri is more close to humans than other species, as this plot might be dominated by strong species specific technical batch effect, rather than biological effect. Moreover, mice and rats seem to be closer to pigs than to themselves from this plot, which is counterintuitive and highlights possible misinterpretation or inadequate analysis.\nResponse: We thank the reviewer for raising this important point. We agree that the MDS panel in Fig. 6C, which separates cells primarily by species, can reflect both true evolutionary divergence and residual technical/batch effects. To mitigate the latter; (1) we limited the comparison to a single tissue (the spleen in all species) to avoid cross-tissue variation; (2) we used a pseudobulk approach (aggregating single-cell data per condition)4 to reduce noise and focus on robust signals; and (3) we restricted the gene set to one-to-one orthologs, ensuring that expression values are comparable across species. Moreover, to minimize technical/“species‑batch” effects, the cross‑species datasets were integrated with Seurat’s reciprocal PCA (RPCA) anchor workflow5, which is designed to conservatively align shared cell states across datasets. After RPCA, UMAP overlays by species show good mixing within matched immune lineages, indicating no obvious technical or species‑driven bias in the integrated space (RFigure 7).\nRFigure 7: Cross‑species spleen immune‑cell landscape. t‑SNE maps of spleen single‑cell transcriptomes from seven vertebrate species.\nImportantly, to avoid heterogeneity across lineages and to align the analysis with our macrophage‑centric conclusions, we now focus the MDS on macrophages only. The revised macrophage‑only MDS provides a more interpretable view that is consistent with the integrated UMAP and our downstream conservation/module analyses (RFigure 8). That said, because MDS is a low‑dimensional projection and cross‑species differences may still capture both biology and residual technical factors, we now treat this panel as qualitative and refrain from over‑interpreting fine‑grained inter‑species distances. Our cross‑species conclusions rely primarily on ortholog‑based comparisons and macrophage‑specific conservation analyses rather than on MDS alone.\nRFigure 8: Multidimensional scaling (MDS) of pseudobulk expression profiles derived from splenic macrophages after one-to-one ortholog mapping. Each symbol represents one species (legend at right). Inter-point distances reflect transcriptome-wide differences and are shown as a qualitative visualization of between-species relationships.\nChanges in the manuscript: Updated the Fig 6C, and revised the Line 330: “To compare cross‑species macrophage programs, we generated pseudobulk profiles by aggregating spleen macrophages across all species to reduce technical noise, and applied multidimensional scaling (MDS) to these profiles. The two‑dimensional MDS projection shows species‑level separation that is broadly consistent with phylogeny; however, because MDS is a projection and cross‑species differences may reflect both biology and residual technical factors, we treat this panel as qualitative and restrict inference to macrophages (Fig. 6C).”\nReviewer #2\nRemarks to the Author:\nThis study by Tang and colleagues presents a comprehensive single-cell immune cell atlas of Tupaia belangeri and provides compelling evidence that the tree shrew is a promising model organism for immunological research. By integrating multi-tissue profiling, EBV infection modeling, and cross-species comparisons, the authors effectively demonstrate the translational potential of Tupaia, particularly in contexts where conventional rodent models are limited.\nImportance: The study is meaningful in that it establishes a foundational resource for using Tupaia belangeri as a non-primate animal model. Its evolutionary proximity to humans and permissiveness to human-specific pathogens make it particularly valuable for translational studies, especially in the field of viral immunology.\nResponse: We sincerely thank the reviewer for the thoughtful and supportive assessment of our study and for highlighting its translational value. We have revised the manuscript to further strengthen the evidence base and to improve clarity and transparency.\nSpecific comments:\n1. Data processing rigor:\n(a) Several aspects of data processing and integration require clearer explanation. For instance, how doublets were removed, what specific QC metrics were applied, and how batch effects were corrected are not sufficiently described. Given the inclusion of multiple tissues and time points, this methodological clarity is essential.\nResponse: We appreciate the reviewer’s request for greater methodological clarity. We have revised the Methods to describe our preprocessing and integration workflow step‑by‑step and to ensure full reproducibility.\n1. Doublet removal. For each library, putative multiplets were identified with DoubletFinder (artificial-nearest-neighbor strategy, default settings) and removed prior to any integration.\n2. Quality control. Cells/nuclei were retained if they had ≥500 detected genes and ≥1,000 UMIs; genes were kept if expressed in ≥5 cells. Cells/nuclei with >30% mitochondrial reads were excluded. To remove outliers, we also discarded the bottom and top 1% of cells by gene and UMI counts per library. Normalization used NormalizeData and the top 2,000 variable genes were selected with FindVariableFeatures (vst).\n3. Batch correction / integration. For the cross-organ atlas (scRNA-seq/snRNA-seq), datasets were integrated using Seurat’s RPCA anchor workflow (FindIntegrationAnchors/IntegrateData). We computed 20 PCs, embedded with UMAP, and clustered with an SNN-Louvain graph (FindClusters, resolution = 1.2). For the EBV spleen time-course (same tissue across conditions), we used SCTransform-based integration (SCT) following Seurat’s recommended pipeline to minimize condition-specific effects while preserving shared cell states.\n4. Annotation. We implemented a two stage procedure: (i) We first separated major classes using canonical markers (e.g., immune vs. non immune; endothelial, epithelial, fibroblast, parenchymal) at the cluster level, guided by cluster defining genes identified with Seurat’s FindAllMarkers (default parameters). This step leverages the marker framework already described in the manuscript and summarized across the 39 cell types. (ii) Immune compartment: we applied CellTypist (model “Immune_All_High.pkl”, default settings) after mapping tree shrew genes to human orthologs, and adopted cluster‑level majority voting with a label‑probability threshold to retain only high‑confidence labels. Non‑immune compartment: we used curated organ‑specific canonical markers from the literature to assign labels and cross‑checked them with the FindAllMarkers outputs. Ambiguous clusters were flagged for manual curation. The final set of 39 annotated cell types and representative marker patterns is summarized in Supplementary Fig. 3, which supports the validity of our labels.\nTo facilitate independent replication, the full analysis pipeline (raw-data processing, QC/doublets, integration, clustering/annotation, cross-species comparisons, figure generation) is publicly available at our GitHub repository (https://github.com/xialabtools/Tree_shrew_scRNA_code). This repository has been archived at Zenodo (https://doi.org/10.5281/zenodo.17609696).\nChanges in the manuscript: Line 522: “Potential doublets were detected and filtered using DoubletFinder63 with default parameters on the basis of the proximity of each cell to artificial doublets. The cells or nuclei were filtered on the basis of a minimal expression threshold of 500 genes and 1000 UMIs, as well as genes that were expressed in at least five cells or nuclei. To remove outliers, we also discarded cells in the bottom 1% and top 1% of the gene and UMI distributions per library. In addition, cells or nuclei with more than 30% mitochondrial gene counts were removed. Filtered data were normalized with Seurat’s normalizeData function and highly variable features (n = 2000) were selected using the “vst” method in “FindVariableFeatures”.Datasets across organs and modalities (scRNA‑seq/snRNA‑seq) were integrated using Seurat’s “the robust principal component analysis (RPCA)” anchor‑based workflow (20 PCs, UMAP embedding; SNN-Louvain clustering at resolution 1.2). For the EBV spleen time‑course we used SCTransform‑based integration (SCT). Cell types were annotated in two stages. First, major cell classes (immune, endothelial, epithelial, fibroblast, and parenchymal compartments) were assigned at the cluster level on the basis of canonical markers, informed by cluster‑defining genes identified with Seurat’s FindAllMarkers (default parameters). Second, immune‑cell subtypes were assigned using CellTypist (model: “Immune_All_High.pkl”; default settings) after mapping T. belangeri genes to their human orthologs. Final labels were derived by a cluster‑level majority‑voting rule subject to a pre‑specified label‑probability threshold. Cells not meeting this threshold were retained as “unassigned” and subsequently curated on the basis of canonical markers. Non‑immune subtypes were annotated using curated organ‑specific marker panels from the literature and cross‑validated against the FindAllMarkers outputs. Representative marker profiles for the 39 annotated cell types are shown in Supplementary Fig. 3.\n(b) More specifically, the localization of T cells in Supplementary Fig. 1 appears questionable, suggesting that doublet removal may not have been adequately performed. Epithelial cell annotations also appear inconsistent across figures. I recommend the authors include UMAP plots annotated by major cell types (e.g., T cells, B cells, epithelial cells, fibroblasts) to improve transparency.\nResponse: We appreciate the reviewer’s careful reading and agree that transparent annotation of major lineages is essential. To address these points, we implemented the following clarifications and additions. Using Seurat’s AddModuleScore (default settings), we computed per‑cell lineage scores for canonical marker panels to visualize and validate major cell types. The marker panels were(RFigure 9):\nB : CD79A, CD19, MS4A1, JCHAIN, XBP1\nT: CD3D, CD3E, CD3G, CD4, CD8A, TYMS\nMyeloid: FCER1G, CSF1R, CD68, KLF1, HNRNPAB, RHCE, ALAS2, GYPA, FCER1A, KIT, PGLYRP1, S100A8, NKG7, LAMP3, CD63\nEndothelial: EGFL7, CEMIP2, FLT1, ADGRF5, DDR2, VWF\nFibroblasts: DCN, FGFR1, NRXN1; HSC: CD34, SLPI, CTSG, MPO\nPlatelets: PPBP, MPIG6B\nParenchyma: LAMA2, MYH7B, AQP5, ASGR1, HNF4A, TFF1, SORBS2, MYL2, AGXT2, MUC5AC, NUPR1\nThis analysis showed that the subset within cluster 18 originally colored as T cells had high B‑cell module scores, while lacking T‑cell markers, indicating B‑cell identity rather than doublets (doublets had already been removed with DoubletFinder). We corrected these labels, which affects only a minor fraction and does not alter any downstream results or conclusions. Beyond this fix, major lineages are consistent with their canonical markers across the atlas.\nRFigure 9: Major‑lineage module scores on the integrated UMAP. Each panel shows the same UMAP of all cells; cells not matching the indicated lineage are plotted in light gray. Color indicates the per‑cell module score computed with AddModuleScore for the specified gene set (higher score = deeper red). Axes: UMAP_1 and UMAP_2; scale bars denote relative module score (arbitrary units).\nTo address the second concern: annotation consistency for epithelial cells and overall transparency. We performed independent Seurat clustering and manual annotation for each organ. This minimizes cross‑tissue transcriptional heterogeneity and refines lineage calls in the appropriate biological background. The per‑tissue UMAPs and marker plot are provided in Supplementary Figs. 19. In addition, We added a Sankey diagram that links L1 major classes → L2 lineage‑level subtypes → L3 tissue‑specific parenchymal identities (e.g., hepatocytes, surface mucous epithelium, proximal tubule, ventricular cardiomyocytes). The clear, predominantly one‑to‑one flows from L1 to L3 demonstrate that organ‑specific parenchymal populations map coherently to their parent lineages, supporting the internal consistency of our annotation strategy (RFigure 10). Finally, we now provide an interactive browser allows readers to query any tissue/cell type/gene for further verification (https://stmm-lab.shinyapps.io/tbca/).\nRFigure 10: Hierarchical mapping of parenchymal identities across the atlas. Sankey diagram linking the L1 major class (Parenchyma; left) to L2 lineage‑level subtypes (middle; Ciliated, Hepatocytes, Surface mucous, Proximal tubule, Ventricular cardiomyocytes) and then to L3 tissue‑specific identities (right). Ribbon width reflects the relative number of cells connecting categories. The clear, predominantly one‑to‑one flows from L1 to L3 show that organ‑specific parenchymal populations map coherently to their parent lineages, supporting the internal consistency of the annotation strategy described in the Methods. Only parenchymal tracks are shown; non‑parenchymal lineages are omitted for clarity.\nChanges in the manuscript: Corrected the small subset in cluster 18 from T to B (Fig. 1C and Supplementary Fig. 1A). Added major-cell-type UMAPs with marker gene set score(Supplementary Fig. 1B ). Added per-tissue UMAPs and annotations (Supplementary Figs. 19).\n(c) Furthermore, batch correction and integration should be benchmarked using multiple methods (e.g., Harmony, scVI) rather than relying on a single approach.\nResponse: We appreciate this suggestion and agree that integration robustness should be demonstrated across methods. In addition to the RPCA‑anchor workflow used in the main analysis, we re‑integrated the dataset with Harmony and FastMNN under matched preprocessing (same HVGs, PCs and clustering parameters). Across the three approaches we observed nearly identical global structure and cell‑type separation, with no spurious segregation by tissue/timepoint and high concordance of lineage markers; only a very small fraction of boundary cells showed method‑sensitive assignment. These results support that our conclusions do not depend on a single integrator. We also considered scVI as requested. Given our data scale (~187k cells across 12 tissues) and our CPU‑only compute environment, training times were prohibitive, so we did not adopt scVI for the production atlas. For the reviewer’s reference, we have included the UMAP plots from the Harmony and FastMNN analyses in this response letter to illustrate the consistency across methods (RFigure 11 and RFigure 12). Since these benchmarks simply reaffirm our original results without altering any conclusions, we have not made changes to the main manuscript, which continues to use the RPCA-based integration pipeline as the core method.\nRFigure 11: Harmony benchmarking of integration. UMAPs of the integrated atlas generated with Harmony under a matched preprocessing setup (doublet removal per library, identical QC thresholds, 2,000 HVGs, 20 PCs) used for the RPCA pipeline. Panels show (top left) Tissue: cells colored by organ of origin; (top right) Cell type: fine‑grained annotations; (bottom left) Major cell type: collapsed lineage‑level classes; and (bottom right) harmony_clusters: SNN–Louvain clusters computed on the Harmony‑integrated embeddings. Batches are well mixed and major biological structures are preserved, yielding cluster and annotation patterns concordant with the RPCA results.\nRFigure 12: FastMNN benchmarking of integration. UMAPs of the integrated atlas generated with FastMNN using the same preprocessing and parameters (doublet removal, QC thresholds, 2,000 HVGs, 20 PCs) as above. Panels show (top left) Tissue, (top right) Cell type, (bottom left) Major cell type, and (bottom right) mnn_clusters (SNN–Louvain on FastMNN embeddings). As with Harmony and RPCA, cell populations remain clearly resolved with consistent cluster structure and biologically coherent global organization, supporting the robustness of our integration strategy.\n2. Depth of analysis in the EBV model:\n(a) The EBV infection model is an important part of the study, but the analysis feels relatively shallow. For example, was EBV transcript detected in B cells or other cell types? Basic analyses such as pseudotime trajectory, cell type-specific differential gene expression (DEG), and pathway enrichment should be provided before diving into cell–cell communication analysis.\nResponse: We thank the reviewer for this insightful comment and would like to clarify that our EBV challenge was designed primarily to induce an in vivo inflammatory environment for studying NR1H3⁺ macrophages, rather than to model EBV pathogenesis itself. In fact, our prior work established that tree shrews can be infected with EBV and mirror key aspects of human EBV infection (e.g., transient acute viremia followed by latent viral persistence)6,7. We built the present study on this established tree shrew EBV model to stimulate the immune system and highlight macrophage dynamics in an EBV-induced inflammatory context.\nTo strengthen the analysis, we have increased the number of animals per time point (now n = 3, 4, and 4 at baseline and two post-infection time points, respectively) to improve statistical power. We also added new in vivo and in vitro validation experiments (see updated Fig. 6 and Supp. Fig. 18) that confirm the activation and characteristics of the NR1H3⁺ macrophages in the EBV-challenged condition, bolstering the conclusions of the study.\nWe appreciate the reviewer’s suggestions for deeper analysis of the single-cell EBV infection data. In response, we performed additional cell-type-specific differential gene expression (DEG) and pathway enrichment analyses with a focus on T cell populations (RFigure 13). These new analyses revealed a robust IFN-γ (IFNG) induction in T cells and an enrichment of a CD8⁺ T cell subset expressing CXCL8, indicative of an active antiviral T cell response. We also observed significant upregulation of viral response pathways at 3 days post-infection (dpi), consistent with early immune activation against EBV. These findings provide further context for the cell–cell communication results and underscore that T cells were engaged during early EBV infection.\nRegarding the detection of EBV transcripts, we re-aligned our scRNA-seq reads to a combined Tupaia belangeri (tree shrew) + EBV reference genome (Cell Ranger pipeline). No EBV transcripts were detected in any cell type, even in B cells. We believe this lack of viral RNA detection is due to several factors, including: (1) Early time points: Our analysis focused on early infection (e.g., 3 dpi), which may precede robust lytic viral gene expression or abundant latent transcripts. (2) Tissue context: The spleen may not be the primary site of EBV replication at these early time points. Given our infection route (intravenous inoculation), the virus might primarily infect peripheral B cells or other tissues before seeding the spleen, resulting in low viral load in the sampled organ. (3) Technical constraints: Droplet-based scRNA-seq has limited sensitivity for low-abundance transcripts, especially for a large DNA virus like EBV. Many EBV transcripts (particularly latent RNAs like EBERs) are not polyadenylated or are expressed at low levels, making them difficult to capture with standard scRNA-seq library preparation and sequencing depth.\nFinally, while we agree that analyses such as pseudotime trajectory mapping of infection could yield interesting insights, we respectfully note that a full trajectory analysis of EBV infection is beyond the scope of our macrophage-centered study. Our focus was to identify and characterize the NR1H3⁺ macrophages under an inflammatory condition induced by EBV, rather than to exhaustively chart the course of EBV infection per se. We consider the reviewer’s suggested trajectory and cell-type resolved analyses very valuable for understanding EBV pathogenesis and immune dynamics; and these will be better addressed in a dedicated future study.\nRFigure 13: Activation and gene expression dynamics of CD8⁺ T cell subsets during early EBV infection. (A) UMAP visualization of T cell subclusters colored by top marker (e.g., CXCL8⁺, CTLA4⁺, GZMH⁺ CD8⁺ T cells). (B) Dot plot showing average expression (color) and percent expression (dot size) of canonical immune checkpoint and cytotoxicity genes across T cell subsets. (C) Barplot showing the relative abundance of selected T cell subsets at 0, 3, and 7 dpi. (D) Differential abundance (cell proportion) estimated with miloR: for each subcluster, neighborhood‑level changes in cell abundance are summarized as log₂ fold change versus 0 dpi; points denote miloR neighborhoods. (E) UMAP overlays of IFNG, TNF, CXCL8, and GZMH expression at each time point (0, 3, 7 dpi), highlighting dynamic cytokine and effector molecule induction. (F) KEGG pathway enrichment (ORA) of DEGs in CD8⁺CXCL8⁺ T cells at 3 dpi, showing significant enrichment of viral response and immune activation pathways.\n(b) Additionally, recent tools like Milo could be combined with improved integration methods to track graph-level changes in cell proportions.\nResponse: We agree and have implemented this analysis. First, the revision increases biological replicates (0/3/7 dpi = 3/4/4), making proportion estimates comparable across time. We then applied miloR to the kNN graph built to test graph‑level differential abundance (DA) of neighbourhoods (3 dpi vs 0 dpi). Milo recapitulates the trends seen in cell proportional changes and immunofluorescence (IF) counts (RFigure 14). Accordingly, in the main text we report the IF quantification as the primary result, because IF directly enumerates NR1H3⁺ and IRF8⁺ macrophages in situ (absolute cells per area) and is agnostic to dissociation and integration choices. The miloR analysis is provided as a robustness check in the Supplement and leads to the same biological conclusion.\nRFigure 14: Cell abundance changes during EBV challenge. (A) Stacked bar plots showing relative proportions of major immune lineages in spleen at 0, 3, and 7 dpi (cohort 3/4/4). (B) miloR neighbourhood graph constructed on the integrated kNN space; each node is a neighbourhood (size ∝ cells, grey edges depict overlap). Node color encodes log₂ fold‑change in neighbourhood abundance (3 dpi vs 0 dpi). (C) Lineage‑wise summary of miloR differential‑abundance statistics. (D) Immunofluorescence (IF) validation: representative spleen sections stained for CSF1R (green) with NR1H3 (red) or IRF8 (yellow) at 0, 3, and 7 dpi; right, absolute counts (cells·mm⁻²) of NR1H3⁺ and IRF8⁺ macrophages per group (mean ± s.e.m.; ns, not significant where indicated).\n(c) Minor comment: In Supplementary Fig. 13A, the temporal order should be displayed as Control → Day 3 → Day 7. Also, the terms “Pos-macro” and “Neg-macro” are not clearly defined and are inconsistent with previous labeling.\nResponse: Thank you for pointing this out. We have corrected Supplementary Fig. 13A to display the time course in chronological order (Control → Day 3 → Day 7). We also replaced the ambiguous labels “Pos-macro” and “Neg-macro” with the consistent terminology used elsewhere in the manuscript: IRF8⁺ macrophages (IRF8_Macro) and NR1H3⁺ macrophages (NR1H3_Macro). These edits standardize nomenclature with the Results section describing NR1H3_Macro and IRF8_Macro (see Fig. 4 text and references to Supplementary Fig. S13 in the manuscript).\nChanges in the manuscript: Supplementary. Fig. 13A panels reordered to Control → 3 dpi → 7 dpi; legend updated accordingly. “Pos‑macro”/“Neg‑macro” replaced with IRF8⁺ macrophages (IRF8_Macro) and NR1H3⁺ macrophages (NR1H3_Macro) in Supplementary Fig. 13.\n3. Cross-species comparison methodology:\n(a) The cross-species integration is one of the more exciting aspects of this paper, but its current implementation is too preliminary. Key technical steps—such as ortholog selection, normalization across species, and gene filtering—should be explained in more detail. Methods like SATURN or scGLUE could be applied here to improve both the interpretability and robustness of the analysis.\nResponse: We thank the reviewer for emphasizing the importance of methodological transparency in the cross‑species analysis. We have expanded the Methods to detail ortholog selection, normalization across species, and the integration workflow, and we provide all code and ortholog tables for full reproducibility.\n1. Focused spleen tissue: To ensure functional comparability, we limited the cross-species analysis to spleen only, thereby reducing tissue heterogeneity across species.\n2. Ortholog selection (gene filtering): Integration was restricted to genes with one-to-one orthologs shared by all species, identified via reciprocal best-hit BLASTP search. This filtering guarantees that only directly comparable genes were used across the datasets. This workflow is consistent with the practice for cross-species single-cell comparisons and has been used in similar contexts (e.g., Tarashansky et al., 20212). A complete ortholog table is provided in the Supplementary Table7.\n3. Normalization and integration: We applied Seurat’s SCTransform (SCT) normalization to each dataset, then performed integration using Seurat’s standard SCT + canonical correlation analysis (CCA) pipeline: SelectIntegrationFeatures → PrepSCTIntegration → FindIntegrationAnchors (with normalization.method = “SCT”) → IntegrateData. The integrated UMAP visualization shows that cells cluster by celltype rather than by species, indicating batch effects arising from species differences were effectively corrected. In other words, different species cells overlap in the integrated embedding, confirming that the integration preserved biological signal (cell identity) while correcting technical differences between species (RFigure 15). The specific script have now released a public GitHub repository (https://github.com/xialabtools/Tree_shrew_scRNA_code).\n4. Method validation: Anchor‑based SCT+CCA has been widely used in high‑impact cross‑species/cross‑cohort atlases and provides a robust solution for single‑modality RNA data, as in our study. Similar workflows have been used in high-impact single-cell studies, including a cross-primate cell atlas8 (Jung et al., Nature 2022) and a human immune cell atlas1 (Domínguez Conde et al., Science 2022), as well as in cross-species tumor microenvironment analyses9 (Elyada et al., Cancer Discovery 2019). If the reviewer or readers feel that employing the newer SATURN tool is still necessary, we are open to exploring it in future work. However, given the scope of our current study (focusing on one organ across species) and the successful outcome with CCA, we believe our chosen approach is sufficient for this analysis.\nRFigure 15: Cross‑species spleen immune‑cell landscape. t‑SNE maps of spleen single‑cell transcriptomes from seven vertebrate species.\nChanges in the manuscript: In the Methods section, we have added a new paragraph under the heading “Cross-species comparisons” (Line: 609):\n“For cross‑species analysis we restricted integration to spleen to minimize tissue‑driven heterogeneity and maintain consistency with the rest of the study. Orthologs were defined by reciprocal best‑hit BLAST of translated peptides against the human proteome; only one‑to‑one orthologs shared by all species were retained and gene identifiers were harmonized to human symbols (Supplementary Table7). For each species, count matrices were processed in Seurat with SCTransform normalization, then performed integration using Seurat’s canonical correlation analysis (CCA) with defualt parameters. For between‑species transcriptomic similarity, we computed per‑species macrophage pseudobulk profiles (log‑normalized within the ortholog set, gene‑wise z‑scored), calculated Euclidean distances, and applied classical multidimensional scaling (MDS) to obtain the first two dimensions shown in the MDS plot.”\n(b) In particular, the MDS plot is difficult to interpret, and it is unclear what the axes represent. If this plot is intended to convey evolutionary distance or transcriptional similarity, it should be supplemented with more formal statistical support.\nResponse: We thank the reviewer for flagging this ambiguity. The classical multidimensional scaling (MDS) panel is intended as a qualitative ordination of transcriptional dissimilarity, not as an evolutionary distance. In other words, they depict relative transcriptional similarities among the species’ macrophage profiles – not a direct evolutionary distance metric. Each point represents the transcriptional profile of macrophages from one species, derived from pseudobulk gene expression aggregates (following the approach of Crowell et al., 20204), thereby minimizing the impact of technical noise and cellular outliers. Then, pairwise Euclidean distances (on log‑normalized, gene‑wise z‑scored expression) were then embedded by classical MDS. The axes (MDS1 and MDS2) in MDS plot are the first two dimensions from classical multidimensional scaling of inter-sample expression distances, analogous to principal components, which are data‑derived coordinates with no intrinsic units. Because MDS is a projection and cross‑species differences may reflect both biology and residual technical factors, we interpret this panel cautiously and restrict inference to spleen macrophages. In additionlly, we have refocused the MDS analysis on macrophages only, rather than displaying all four immune lineages together. This targeted approach reduces heterogeneity from other lineages and aligns the figure with our macrophage-centered narrative.\nChanges in the manuscript:\n1. Fig. 6C legend updated: “Multidimensional scaling (MDS) plot of pseudobulk transcriptomes showing the relationships between different species. Each point denotes the per‑species spleen macrophage expression profile. Pairwise Euclidean distances were computed on log‑normalized, gene‑wise z‑scored expression restricted to one‑to‑one orthologs shared across species; coordinates (MDS1/MDS2) are from classical MDS (cmdscale). Axes visualize relative transcriptional similarity and have no intrinsic units; the panel is intended as a qualitative ordination rather than an evolutionary distance”.\n2. Methods (Cross‑species comparisons) now includes the exact MDS computation (Line 577): “For between‑species transcriptomic similarity, we computed per‑species macrophage pseudobulk profiles (log‑normalized within the ortholog set, gene‑wise z‑scored), calculated Euclidean distances, and applied classical multidimensional scaling (MDS) to obtain the first two dimensions shown in the MDS plot”.\n3. Results interpret cautiously (Line 330): “To compare cross‑species macrophage programs, we generated pseudobulk profiles by aggregating spleen macrophages across all species to reduce technical noise, and applied multidimensional scaling (MDS) to these profiles. The two‑dimensional MDS projection shows species‑level separation that is broadly consistent with phylogeny; however, because MDS is a projection and cross‑species differences may reflect both biology and residual technical factors, we treat this panel as qualitative and restrict inference to macrophages (Fig. 6C).”\n(c) Additionally, the analysis of HLA loci, which are known to be highly divergent between species, needs to be handled with extra caution. Ortholog mapping for HLA genes should be genome-informed and well-documented.\nResponse: We thank the reviewer for highlighting the special challenges of cross‑species analysis at the HLA (MHC) loci. We fully agree that ortholog assignment in this region is complicated and require special caution by lineage‑specific duplications/deletions and incomplete one‑to‑one correspondence across species10. To avoid over‑interpretation, we have removed all cross‑species, gene‑level analyses involving classical MHC/HLA loci from the manuscript. In the revised cross‑species comparisons, we have omitted the HLA ortholog comparison, focusing instead on more conserved immune genes and pathways where orthology is clear and well-supported. Analyses that involve HLA expression are retained only within human‑specific contexts and are no longer used for cross‑species inference.\nChanges in the manuscript: The cross-species HLA analysis has been removed.\nReviewer #3\nRemarks to the Author:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We thank you for co‑reviewing our manuscript as part of Nature Communications’ Early Career Researcher peer‑review initiative. We are grateful for the time and care you and your co‑reviewer devoted to this review; your constructive comments materially improved the clarity and rigor of our work. We fully support this initiative and appreciate the opportunity to both benefit from and contribute to ECR training in peer review.\nReviewer #4\nRemarks to the Author:\nThis new study describes a single cell immune atlas for the tree shrew (Tupaia belangeri). Tree shrews are an emerging small animal model that are evolutionarily closer to primates than mice. As such, tree shrews represent a useful, viable alternative model for infectious disease and aging studies that may be otherwise challenging to carryout in murine or non-human primate model systems. The investigators profiled single cells from 12 different tissues of 3 adult tree shrews, identifying 39 different cell types. Notably, bulk RNAseq studies have been performed in spleen, PBMCs, and other tissues from the tree shrew, and single cell transcriptome atlases have been recently published for brain and neurological tissues. A strength of this study is that it represents the first tree shrew immune-focused single cell atlas which can provide an important platform for understanding systemic transcriptional responses to infectious diseases in future studies. Aside from the substantial efforts that went into generating this resource, there are some key weaknesses and major concerns with the overall conclusions derived from the interpretation of the EBV infection data.\nResponse: We thank Reviewer for recognizing the value of the first immune‑focused single‑cell atlas of Tupaia belangeri and for the constructive critique of our EBV analysis. We have revised manuscript according to your suggestion. Our explanation to the comments point- by- point is as follow.\nMajor Concerns:\n1. Atlas reveals that T.belangeri possess well-developed tonsillar structures analogous to human tonsils. These tonsillar anatomical structures within the oropharynx were previously reported for the tree shrew—and thus, it is expected that these tissues will harbor lymphoid cell types commonly observed in primate tonsils.\n(a) In figure 1, mature and activated B cell subsets are not detected in the pharyngeal tonsils, which is unusual. Could this be due to sampling issues or issues with adequate sequencing coverage?\nResponse: We agree that the apparent absence of mature and activated B cell subsets in the pharyngeal tonsil is unexpected. We believe this reflects sampling limitations rather than sequencing coverage. The pharyngeal tonsil is located deep near the skull base (RFigure 16) and was technically difficult to harvest; despite intending to sample three animals, only one specimen yielded sufficient tissue for scRNA‑seq. With this limited input, rare lymphocyte states are likely under‑represented even though per‑cell quality was adequate.\nNotably, we still detected B cells (6.7%) and plasma cells (3.7%) in the pharyngeal tonsil. To verify that germinal‑center B cells were present in the tissue but missed at low frequency in scRNA‑seq, we performed immunohistochemistry, which showed BCL6⁺ (germinal‑center) and CR2/CD21⁺ B cells in pharyngeal‑tonsil sections (Supplementary Fig. 11). We have clarified in the text that the tonsil data represent a baseline anatomical/cellular atlas; the lower representation of mature/activated B‑cell subsets in this particular pharyngeal sample is most plausibly a sampling issue rather than a biological absence.\nRFigure 16: Surgical collection of tree-shrew pharyngeal tonsil. (A–D) Step-wise exposure of the nasopharyngeal cavity after median mandibular split. (E) Endoscopic view of the pharyngeal tonsil (dashed outline) on the posterior nasopharyngeal wall before excision. (F) Resected specimens: lower fragment, pharyngeal tonsil; upper fragment, soft palate (scale in cm).\nChanges in the manuscript: Results (Anatomical and functional validation of tonsillar immunity in T. belangeri) now includes the clarifications (Line 184): “In the pharyngeal tonsil, scRNA‑seq captured only small numbers of B cells and plasma cells and did not recover mature/activated B‑cell subsets; however, immunohistochemical analysis demonstrated germinal centers, indicating these populations exist in vivo and were likely missed due to limited tissue yield rather than inadequate sequencing coverage.”\n(b) As the tonsils are a significant focus of the paper (Fig 2), more work needs to be put into describing these cells. Fig. 2B needs to be described in detail—UMAP represents cells in which tonsils?\nResponse: We thank the reviewer for this helpful suggestion and agree that additional clarity was needed. Figure 2B depicts the UMAP of single cells from the palatine tonsil only; the pharyngeal tonsil was analyzed separately and is presented in Supplementary Fig. 11. We have revised the main text and figure legend accordingly and expanded the description of Fig. 2B to specify the principal clusters and example markers, including B‑cell lineages (CD79A/MS4A1; plasma cells MZB1/JCHAIN), CD4⁺ T‑cell states (naïve LEF1/SELL, memory S100A4/TNFRSF18, Treg FOXP3/LAG3, exhausted CTLA4/TOX), myeloid populations (macrophages CSF1R/CD68; dendritic cells FCER1A/ITGAX), and epithelial/secretory cells (EPCAM/KRT14/KRT16).\nChanges in the manuscript: Results (Anatomical and functional validation of tonsillar immunity in T. belangeri) now includes the sentence (Line 147): “Single-cell sequencing of the palatine tonsils revealed 17 distinct cell populations span immune and epithelial/stromal compartments, including B cells (CD79A/MS4A1; plasma cells MZB1/JCHAIN), CD4⁺ T‑cell states (naïve LEF1/SELL, memory S100A4/TNFRSF18, Treg FOXP3/LAG3, exhausted CTLA4/TOX), myeloid populations (macrophages CSF1R/CD68, dendritic cells FCER1A/ITGAX), and epithelial/secretory cells (EPCAM/KRT14/KRT16) (Fig. 2B-C). Pharyngeal tonsil cells were analyzed separately and are shown in Supplementary Fig. 11.”\nFigure 2 legend (first line): “Fig. 2B: UMAP of palatine tonsil single cells; the pharyngeal tonsil is analyzed separately (Supplementary Fig. 11).”\n(c) There were only ~2700 cells sequenced for pharyngeal-were palantine and pharyngeal grouped together in this figure?\nResponse: We appreciate the reviewer’s question. No—palatine and pharyngeal tonsil datasets were analyzed and visualized independently throughout. As shown in Fig. 1B, the pharyngeal tonsil dataset contained ~2,700 cells, whereas the palatine tonsil dataset comprised ~14,700 cells; accordingly, Fig. 2B displays only the palatine tonsil UMAP, and the pharyngeal tonsil results are presented in Supplementary Fig. 11. We have clarified this in the text and figure legend.\nChanges in the manuscript: The pharyngeal tonsil yielded approximately 2,700 cells (as reported in Line 98 of the manuscript), whereas the palatine tonsil yielded a much larger number of cells (14701 cells). We confirm that we analyzed and visualized the palatine and pharyngeal tonsil datasets separately at all times. Figure 2B represents the palatine tonsil cells only. Pharyngeal tonsil cells were analyzed separately and are shown in Supplementary Fig. 11.\n(d) Are the cell subsets and relative distributions of representative cell types in these tissues comparable to what is reported in human and nonhuman primate samples?\nResponse: We agree that it is important to benchmark tree‑shrew tonsillar biology against human/NHP data. Across three orthogonal evidence, we find broad concordance:\n1. Anatomy and H&E morphology. The tree shrew (T. belangeri) palatine tonsils are located in the lateral oropharynx (RFigure 17), and a pharyngeal tonsil resides in the nasopharyngeal roof, forming a Waldeyer’s ring equivalent, just as in primates. These tonsillar tissues lie adjacent to the airway and digestive tract, analogous to human tonsillar anatomy and differs from rodents like rodents that lack true tonsils11. Histological sections reveal numerous lymphoid follicles with germinal centers (GCs) beneath specialized epithelia, mirroring the fundamental organization in human tonsils12. In the palatine tonsil, a stratified squamous epithelium (non-keratinized) covers deep invaginating crypts, underlain by dense follicles. This is identical to human palatine tonsils, which feature crypt epithelium and follicle-rich lymphoid tissue for antigen trapping13. We also note the germinal center morphology in tree shrew tonsils is highly characteristic: pale GC zones with dark peripheral mantles can be seen (H&E), indicating active B-cell proliferation and selection, exactly as in human. The overall compartmentalization with B cell follicular zones and T cell-rich interfollicular areas is conserved. In summary, the tree shrew tonsils exhibit the same tissue organization, epithelial characteristics, and lymphoid microarchitecture described in human and NHP tonsils.\nRFigure 17: Histology of tree‑shrew tonsils (H&E). (A, H) Low-magnification transverse sections through the naso-/oropharyngeal region showing the pharyngeal tonsil in the nasopharyngeal roof and palatine tonsils in the lateral oropharynx; boxed areas indicate regions shown at higher magnification. (B–G) Higher-magnification views from panel A highlighting lymphoid follicles with germinal centers beneath specialized epithelia: stratified squamous crypt epithelium and ciliated respiratory epithelium overlying follicle-rich lymphoid tissue. (I–N) Higher-magnification views from panel H demonstrating additional micro-architectural features, including interfollicular zones, crypt structures, and epithelial specializations (e.g., ciliated mucosa and non-keratinized squamous epithelium) adjacent to follicles.\n2. Protein‑level zoning (IHC). Marker localization recapitulates the canonical tonsillar micro‑architecture (RFigure 18). For example, T cells marked by CD3D are localized to the interfollicular T-cell zones and around crypts, whereas B cells (marked by CD79A, MS4A1/CD20, etc.) densely populate the follicles. This compartmentalized pattern of CD3+ T-cell zones adjacent to CD79A+ B-cell follicles is a hallmark of human tonsil histology. Furthermore, within germinal centers of tree shrew, we find the expected proliferation and germinal-center markers: many GC B cells are strongly MKI67 (Ki-67) positive, and express the transcription factor BCL6 in their nuclei, reflecting active centroblast proliferation and GC identity. The prevalence of Ki-67 high cells in tree shrew GCs (~80% of GC B cells) is on par with that in normal human tonsils. We also examined CR2 (CD21), a complement receptor highly expressed on follicular dendritic cells (FDCs) that form the stromal meshwork in GCs. In tree shrew GCs, CR2+ FDC networks are readily detected, outlining the GC structure – exactly as in human tonsils, where FDCs universally express CD21 in the light and dark zones14. We note that epithelial cell markers (e.g. pan-cytokeratins) in the tree shrew tonsil also showed expected patterns (highlighting the stratified squamous epithelium of palatine tonsil and ciliated respiratory epithelium of pharyngeal tonsil), further underscoring conserved tissue specialization. In summary, the spatial localization of immune cell subsets in the tree shrew tonsil (T cells, B cells, germinal center B cells, and FDCs) recapitulates the known histological organization of human tonsils, including the proper segregation of markers like CD3, CD20/CD79a, Ki-67, BCL6, and CD21.\nRFigure 18: Immunohistochemistry of tree-shrew tonsillar tissues. Representative IHC of the nasopharyngeal (pharyngeal) tonsil (top block) and palatine tonsil (bottom block). Panels show: CD3 (T-cell zones in interfollicular areas/around crypts), CD79a (B-cell follicles), Ki-67 (proliferating cells within germinal centers), BCL-6 (GC B-cell nuclei), CD21 (follicular dendritic-cell networks outlining GCs), and epithelial markers EPHA2 (surface/crypt epithelium) and AQP5 (submucosal gland/duct epithelium). Scale bars as indicated.\n3. scRNA‑seq signatures: Single-cell data also show clear concordance of canonical tonsillar cell markers between tree shrews and humans. In the tree shrew palatine tonsil, we identified the expected immune cell subsets by hallmark genes – B cells (CD79A, MS4A1), plasma cells (JCHAIN, MZB1, XBP1), germinal-center B cells (BCL6, MKI67, AICDA), naïve/memory CD4+ T cells (LEF1, SELL, S100A4, TNFRSF18), Treg cells (FOXP3, LAG3), exhausted CD4+ T subset (CTLA4, TOX), CD8+ T/NK lymphocytes (CD8A, GNLY, NKG7), myeloid cells (CSF1R, CD68), and epithelial/secretory lineages (EPCAM, KRT14, AQP5, EPHA2). The pharyngeal tonsil (adenoid) contained parallel populations, including ciliated respiratory epithelial cells (TPPP3, NEK10), glandular cells (AQP5, FAM3D), goblet-like cells (KRT15, VMO1), and abundant lymphocytes (e.g. CD3D+ T cells, CD79A/CD19+ B cells, and JCHAIN+ plasma cells). These marker expression patterns are concordant with recent human tonsil single‑cell atlases (e.g., Massoni‑Badosa et al., Immunity 2024) that annotate the same compartments using these canonical genes15. We have added a brief statement in the Discussion noting the concordance of tonsillar architecture and canonical cell‑type markers between tree shrew and human tonsils, and we cautiously frame this as evidence of the tree shrew’s potential rather than superiority over traditional rodent models for studying pharyngeal lesions and respiratory‑acquired infections.\nChanges in the manuscript: We have added a sentence description in the Discussion Line 396: “For example, our data show the palatine and pharyngeal tonsils of T. belangeri lie at mucosal interfaces directly exposed to the external environment and therefore constitute a first line of defence.. At the molecular level, T. belangeri tonsils display canonical marker programs across B‑cell (including GC), T‑cell, myeloid, and epithelial/secretory compartments (e.g., CD79A/MS4A1; JCHAIN/MZB1; BCL6/MKI67/AICDA; FOXP3/LAG3; CTLA4/TOX; CD8A/GNLY/NKG7; EPCAM/KRT14; AQP5/EPHA2), consistent with the human tonsil atlas of Massoni‑Badosa et al33. Together with the presence of well‑formed germinal centres and abundant lymphocytes which features not observed in rodent models48. These findings indicate similarity in immune landscape between T. belangeri and humans. Accordingly, T. belangeri shows promising application for studying pharyngeal lesions and respiratory‑acquired infectious diseases.”\n2. A major caveat is that the EBV infection model does not mimic natural oral uptake and natural viral transmission which is observed in other NHP models with LCVs (i.e. rhesus and cynomolgus macaque models). Here, the route of EBV administration is intravenous into the femoral vein, thus any relevance to palantine/pharyngeal tonsil responses—which the authors argue is an important component of the tree shrew model (Fig. 2)—is unclear as this requires trafficking of the virus to the oral/throat cavity (in previous publications with this tree shrew model, EBV shedding occurs at 7 dpi). It is unclear why tonsil cells were not sequenced from the infected animals as this would highlight the importance and biological relevance of these structures for EBV infection. The I.V. infection route and also the high titer virus administered is akin to what is utilized in mouse models, which does not recapitulate natural EBV infection in humans and this needs to be acknowledged. As such, the immune responses (i.e. bypassing mucosal immunity in the oral cavity and gut that occurs with oral challenge) are expected to be skewed.\nResponse: We thank the reviewer for this thoughtful comment and agree that intravenous (i.v.) EBV challenge does not recapitulate natural oral transmission. We have clarified this limitation and our rationale as follows:\n1. Study aim and rationale for i.v. infection. Our goal was to resolve early, spleen‑centered transcriptional dynamics of NR1H3⁺ tissue macrophages. To obtain a synchronous systemic exposure at defined early time points (3–7 dpi), we used a femoral‑vein inoculum. This minimizes the variability in uptake/kinetics inherent to oral challenge and enables detection of transient systemic programs with cellular resolution.\n2. Explicit acknowledgement of route‑related limitations. We fully acknowledge that i.v. administration bypasses the oral mucosa and can skew early immune responses relative to natural infection by minimizing mucosal sensing. We now state this clearly in the Discussion and, in Methods, specify the inoculum strain, dose and route together with the rationale for using i.v. solely to synchronize early systemic exposure. Accordingly, our conclusions are restricted to early systemic (not mucosal) responses.\n3. Why tonsils were profiled only at baseline. We profiled palatine and pharyngeal tonsils from uninfected animals to establish that tree shrews possess organized, human‑analogous tonsillar structures and to provide a reference atlas of their cellular composition (Fig. 2). These baseline datasets support the model’s suitability for future mucosal studies but are independent of the acute i.v.‑challenge cohort.\n4. Why we did not sequence tonsils post‑infection here. The infection arm was designed around early acute time points (≤7 dpi) to interrogate immediate systemic responses in the spleen. Given this design, we did not power the study to capture oropharyngeal trafficking or early tonsillar events; therefore, we cannot draw conclusions about mucosal infection from these data. We now make this explicit in the manuscript.\n5. Future directions: oral challenge and mucosal readouts. In response to the reviewer’s suggestion, we have added to the Discussion a dedicated paragraph outlining plans to establish an oral EBV challenge in tree shrews, with time‑resolved sampling of palatine/pharyngeal tonsils and paired throat‑swab virology, to directly compare mucosal and i.v. routes.\nWe are grateful for this constructive critique; the revision strengthens the framing and clearly delimits the scope of our conclusions.\nChanges in the manuscript: We have added a sentence description in the Discussion Line 481: “Finally, because we used intravenous EBV inoculation to synchronize early systemic responses, our study does not model natural oral transmission or early tonsillar immunity. Mucosal innate and adaptive events may therefore differ from the responses reported here. Future work will establish an oral challenge in tree shrews with time‑resolved tonsillar profiling and shedding/transmission readouts to benchmark against the i.v. paradigm.”\n3. The sample size for EBV experiments in Fig 4 is extremely small. Spleen cells are sequenced from only two animals that were infected with EBV and harvested at different time points, limiting the statistical power and generalizability of findings that NR1H3(LXRa+) macrophages are indeed altered by EBV infection. It is basically a pilot experiment. In an attempt to confirm scRNAseq findings, immunofluorescence and qRT-PCR experiments were performed on samples from 4 animals in Fig 4. No differences in NR1H3+ or IRF8+ macrophage numbers were detected by IFA, although the NR1H3 RNA level does decrease, which hardly confirms the scRNAseq findings. Given that data describing NR1H3 cell reprogramming is from n=1, additional evidence should be provided to support the major conclusions that NRH13+ macrophages dominate the inflammatory responses observed in EBV infection. For example, does depletion of these cells impact acute EBV infection? Are there additional signs of inflammatory markers that track with these macrophage subsets?\nResponse: We understand this concern. We agree that the original EBV arm in which one biological replicate per time point was profiled by scRNA‑seq was underpowered. In the revised study we have taken the following steps to strengthen the evidence and clarify our claims:\n1. Increased animal numbers. We expanded cohort sizes to Control n=3, 3 dpi n=4, 7 dpi n=4 for scRNA‑seq, and to n=6/group for multiplex IFA and qRT‑PCR. The major biological patterns observed in the pilot dataset were reproduced including the early shift of NR1H3⁺ macrophages toward increased network centrality and inflammatory output, thereby addressing concerns about statistical power and generalizability (new Figure 4).\n2. What we mean by “dominant” and why IFA counts do not change. Our conclusion does not hinge on numerical expansion of NR1H3⁺ macrophages. Consistent with the reviewer’s observation, IFA quantification shows stable numbers of NR1H3⁺ and IRF8⁺ macrophages across time points. Instead, the key finding is functional reprogramming: during acute infection, NR1H3⁺ macrophages become major signalling hub, with heightened incoming/outgoing communication strength. This is accompanied by loss of anti‑inflammatory programs in NR1H3_Macro (UCell score) while pro‑inflammatory programs in IRF8_Macro remain unchanged. Transcriptionally, NR1H3 itself decreases and TNF/IL1B/TNFSF and C1Q components increase in NR1H3_Macro. In the revision, we explicitly rephrase “dominates” as “emerges as a principal signalling hub and cytokine source during acute infection,” to avoid implying a change in cell abundance.\n3. Depletion experiments. We agree that subset‑specific depletion in vivo would be definitive. At present, species‑specific depleting antibodies/conditional tools for tree shrew macrophage subsets are not available, precluding clean depletion without broad off‑target effects. We now state this constraint transparently and limit our claim to network‑level dominance and cytokine output rather than causal necessity in vivo. We also outline in the Discussion that future work will pursue targeted perturbations of the NR1H3 pathway once suitable in‑species tools are feasible.\n4. Additional Functional evidence (human macrophage systems and ex vivo validation). To move beyond transcriptional associations, we added gain‑of‑function / pharmacologic tests showing that activating NR1H3 dampens inflammatory output: in THP‑1 macrophages, NR1H3 activation reduces IL‑1β, IL‑6, TNF‑α and diminishes effector polarization/IFN‑γ production of primary CD8⁺ T cells while leaving T‑cell proliferation unchanged; in primary human monocyte‑derived macrophages, the NR1H3 agonist T0901317 similarly suppresses LPS‑induced IL‑1β, TNF‑α, IL‑6 (RFigure 19, now compiled as Supplementary Fig. 18). These results are consistent with our model that down‑tuning NR1H3 removes a brake, enabling pro‑inflammatory output during EBV challenge.\nRFigure 19: Functional of NR1H3 activation in human macrophages and translational ex vivo validation. (A) THP 1–derived macrophages (WT, NR1H3 KO and NR1H3 OE) were stimulated under a standardized pro inflammatory challenge; supernatants were collected and IL 1β, IL 6 and TNF α were quantified by ELISA (Elabscience; OD 450 nm, BioTek). NR1H3 OE released significantly less of all three cytokines than WT and KO, whereas KO showed an increase for IL 1β and no significant change for IL 6/TNF α under these conditions. (B) Phagocytosis assay using pHrodo™ Red E. coli bioparticles (Thermo Fisher). THP 1 cells were differentiated with PMA 100 ng/mL for 24 h, incubated with bioparticles (1:100, 4 h), and analyzed by flow cytometry; NR1H3 OE displayed a reduced pHrodo positive fraction compared with WT. (C) Schematic of the T cell crosstalk workflow: CD3⁺ T cells were negatively selected from healthy donor PBMCs (written consent; Ethics 2025 E0760), PHA stimulated, and cultured with conditioned supernatants from WT or NR1H3 OE THP 1 macrophages. (D) IFN γ ELISPOT showed fewer spots in T cells exposed to NR1H3 OE supernatant than to WT supernatant. (E) Flow cytometry with intracellular cytokine staining demonstrated a lower frequency of IFN γ⁺ cells within CD8⁺ T cells after exposure to NR1H3 OE supernatant. (F) CFSE dilution analysis indicated no difference in T cell proliferation between conditions. (G) Ex vivo validation in primary human monocyte derived macrophages: CD14⁺ monocytes were isolated by negative selection, differentiated for 5 days with M CSF (50 ng/mL) to adherent macrophages, treated with the NR1H3 agonist T0901317 for 24 h, then challenged with LPS as indicated. (H) ELISA of macrophage supernatants showed significant reductions in IL 1β, TNF α and IL 6 after T0901317 treatment compared with control. Data are mean ± SD; each dot represents an independent experiment or donor. *P<0.05; **P<0.01; ns, not significant.\nTogether, the expanded cohorts, clarified wording, and orthogonal functional evidence substantially strengthen the conclusion that NR1H3⁺ macrophages undergo early functional reprogramming and emerge as a principal signalling hub during acute EBV infection, even though their absolute numbers remain stable.\nChanges in the manuscript:\n1. We re‑computed the data integration and downstream analyses on the updated dataset and replaced Figure 4 and its legend accordingly.\n2. The heading was changed from “Dominant role of T. belangeri NR1H3⁺ macrophages in resolving acute viral inflammation” to “Early functional reprogramming of NR1H3⁺ macrophages during acute EBV infection.” Correspondingly, within the Results we replaced “dominate” with “merge as a principal signaling hub/cytokine source,” to avoid implying a change in cell abundance.\n3. We added Supplementary Fig. 18 compiling THP‑1 and primary human macrophage assays (cytokine release, phagocytosis, and T‑cell readouts). In the Results (Line 287), we inserted the following text to reference these data: “In functional assays, NR1H3 activation in THP‑1 macrophages reduced secretion of IL‑1β, IL‑6 and TNF‑α, modestly decreased uptake of bacterial particles, and diminished IFN‑γ production and effector polarization of primary CD8⁺ T cells while leaving T‑cell proliferation unchanged (Supplementary Fig. 18A–F). Extending these findings ex vivo, pharmacologic activation of NR1H3 with T0901317 in primary human monocyte‑derived macrophages similarly suppressed LPS‑induced IL‑1β, TNF‑α and IL‑6 release (Supplementary Fig. 18G,H), underscoring the translational relevance of NR1H3‑mediated immunoregulation.”\n4. Depletion of liver X receptors (NR1H3 gene encodes LXRa) in mice leads to enhanced replication of a murine g-herpesvirus (MHV68) that is often used as a model for EBV infection (PMC6050960). Such prior studies indicate NRH13-mediated mechanisms are indeed conserved in mice.\nResponse: We appreciate the pointer to the MHV68 study showing that loss of NR1H3 (LXRα) augments γ‑herpesvirus replication in mice, supporting a conserved immunoregulatory role for LXRα. Our work builds on and complements these findings in three ways.\n1. We delineate a distinct NR1H3⁺ splenic macrophage subset with robust anti‑inflammatory/immunosuppressive features; this subset has not been explicitly delineated in mouse spleen datasets.\n2. Mechanistically, we show that NR1H3 restrains the non‑canonical NF‑κB program: genetic and pharmacologic perturbations reduce NFKB2, MAP3K14 (NIK), NFKBIZ, and Cut&Tag reveals NR1H3 occupancy at/near regulatory regions of TNFSF14, NFKB2, NFKBIZ.\n3. We provide cross‑species support, identifying a functionally analogous NR1H3⁺ subset in human spleen and showing by comparative single‑cell profiling that this program is readily detected in tree shrew and human, whereas it appears less prominent in mouse spleen, consistent with species‑specific cellular contexts rather than a lack of pathway conservation.\nTogether, these data refine how NR1H3 constrains inflammation via repression of non‑canonical NF‑κB in splenic macrophages and complement mouse studies by revealing a subset‑level, human‑relevant axis that is well captured in the tree‑shrew model.\nChanges in the manuscript: In light of prior mouse work, we softened cross‑species wording throughout and avoided categorical statements (e.g., “absent in mice”) and softened cross‑species phrasing, clarifying that our findings complement the mouse literature.\nLine 35: Deleted sentence: “yet absent in mice.”\nLine 38: Deleted sentence: “that more closely mirrors human immunity than mice”.\nRevised sentence: “Together, these findings highlight a conserved anti‑inflammatory mechanism shared by tree shrew and human, emphasizing the tree shrew’s unique value as an immunological model that complements rodent and primate systems.”\nLine 350: Delete sentence: “and highlighting the potential advantage of using T. belangeri as a more suitable model for human disease research”\nLine 457: Revised sentence: “These findings support T. belangeri as a complementary model for studying immune responses.”\n5. The study largely focuses on describing B cells and macrophage subset, yet T cell responses are a major contributor to control of infectious agents that have been tested in tree shrew models (i.e. Zika, hepatitis virus). T cell populations should be described in detail.\nResponse: We thank the reviewer and agree that T cells are pivotal contributors to antiviral control in tree shrew models. While our study was primarily designed to resolve NR1H3⁺ macrophage programs under acute EBV challenge, we have now added a dedicated, lineage‑resolved analysis of splenic T cells (0/3/7 dpi) placed immediately before the cell–cell communication section (Supplementary Fig. S14). These new analyses revealed a robust IFN-γ (IFNG) induction in T cells and an enrichment of a CD8⁺ T cell subset expressing CXCL8, indicative of an active antiviral T cell response. We also observed significant upregulation of viral response pathways at 3 days post-infection (dpi), consistent with early immune activation against EBV (RFigure 20). These findings provide further context for the cell–cell communication results and underscore that T cells were engaged during early EBV infection. Given our macrophage‑centered objective, we did not pursue full T‑cell trajectory or TCR repertoire analyses here. We agree these will be informative, and we plan a follow‑up, EBV‑focused topic study to probe T‑cell dynamics in depth comparing oral vs i.v. challenge, sampling tonsil and other mucosal sites, and adding clonotype analyses across finer time points.\nRFigure 20: Activation and gene expression dynamics of CD8⁺ T cell subsets during early EBV infection. (A) UMAP visualization of T cell subclusters colored by top marker (e.g., CXCL8⁺, CTLA4⁺, GZMH⁺ CD8⁺ T cells). (B) Dot plot showing average expression (color) and percent expression (dot size) of canonical immune checkpoint and cytotoxicity genes across T cell subsets. (C) Barplot showing the relative abundance of selected T cell subsets at 0, 3, and 7 dpi. (D) Differential abundance (cell proportion) estimated with miloR: for each subcluster, neighborhood‑level changes in cell abundance are summarized as log₂ fold change versus 0 dpi; points denote miloR neighborhoods. (E) UMAP overlays of IFNG, TNF, CXCL8, and GZMH expression at each time point (0, 3, 7 dpi), highlighting dynamic cytokine and effector molecule induction. (F) KEGG pathway enrichment (ORA) of DEGs in CD8⁺CXCL8⁺ T cells at 3 dpi, showing significant enrichment of viral response and immune activation pathways.\nOther minor comments:\n1. Introduction, line 52: smaller NHPs such as marmosets and tamarins have been historically used for EBV infection studies and should be mentioned in this section.\nResponse: Thank you for the suggestion. We have revised the text accordingly.\nChanges in the manuscript: In Introduction Line 52, the following text has been added: “Given their close genetic and immunological proximity to humans, nonhuman primates (NHPs) provide more translatable insights; for example, smaller New World NHPs such as cottontop tamarins and common marmosets have served as experimental systems for EBV infection, including EBV-positive lymphoma models16.”\n2. Line 74, first single-cell atlas of T. belangeri and line 435: This is incorrectly stated as there are at least three other publications describing brain/neurological tissue-focused single cell atlases for the tree shrew.\nResponse: Thank you for pointing this. We have revised the text accordingly.\nChanges in the manuscript: Line 25; Line 395; Line 488: deleted “first”.\nReferences\n1. Domínguez Conde, C. et al. Cross-tissue immune cell analysis reveals tissue-specific features in humans. Science (New York, N.Y.) 376, eabl5197; 10.1126/science.abl5197 (2022).\n2. Stadtmauer, D. J. et al. Cell type and cell signalling innovations underlying mammalian pregnancy. Nature ecology & evolution 9, 1469–1486; 10.1038/s41559-025-02748-x (2025).\n3. Zhang, L. et al. Lineage tracking reveals dynamic relationships of T cells in colorectal cancer. Nature 564, 268–272; 10.1038/s41586-018-0694-x (2018).\n4. Crowell, H. L. et al. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat Commun 11, 6077; 10.1038/s41467-020-19894-4 (2020).\n5. Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573-3587.e29; 10.1016/j.cell.2021.04.048 (2021).\n6. Wang, Z. et al. A study of Epstein-Barr virus infection in the Chinese tree shrew(Tupaia belangeri chinensis). Virology journal 14, 193; 10.1186/s12985-017-0859-5 (2017).\n7. Xia, W. et al. Tree Shrew Is a Suitable Animal Model for the Study of Epstein Barr Virus. Frontiers in Immunology 12, 789604; 10.3389/fimmu.2021.789604 (2021).\n8. Jung, M. et al. Cross-species transcriptomic atlas of dorsal root ganglia reveals species-specific programs for sensory function. Nat Commun 14, 366; 10.1038/s41467-023-36014-0 (2023).\n9. Elyada, E. et al. Cross-Species Single-Cell Analysis of Pancreatic Ductal Adenocarcinoma Reveals Antigen-Presenting Cancer-Associated Fibroblasts. Cancer discovery 9, 1102–1123; 10.1158/2159-8290.CD-19-0094 (2019).\n10. Trowsdale, J. & Knight, J. C. Major histocompatibility complex genomics and human disease. Annual review of genomics and human genetics 14, 301–323; 10.1146/annurev-genom-091212-153455 (2013).\n11. Wu, R.-Q., Zhang, D.-F., Tu, E., Chen, Q.-M. & Chen, W. The mucosal immune system in the oral cavity-an orchestra of T cell diversity. International journal of oral science 6, 125–132; 10.1038/ijos.2014.48 (2014).\n12. Brandtzaeg, P. Immunobiology of the Tonsils and Adenoids. In Mucosal Immunology (Elsevier2015), pp. 1985–2016.\n13. Tonsils and adenoids: everything the entsurgeon needs to know (2002).\n14. Imal, Y. & Yamakawa, M. Morphology, function and pathology of follicular dendritic cells. Pathology international 46, 807–833; 10.1111/j.1440-1827.1996.tb03555.x (1996).\n15. Massoni-Badosa, R. et al. An atlas of cells in the human tonsil. Immunity 57, 379-399.e18; 10.1016/j.immuni.2024.01.006 (2024).\n16. Wang, F. Nonhuman primate models for Epstein-Barr virus infection. Current opinion in virology 3, 233–237; 10.1016/j.coviro.2013.03.003 (2013).\n\nResponse to Reviewers’ Comments\nWe thank the Editor and the reviewers for their careful evaluation of our manuscript and for the constructive comments provided in this round. We have revised the manuscript accordingly and addressed all remaining points raised by the reviewers. In particular, we have clarified the experimental design and data interpretation for the ELISA and CFSE assays, corrected figure labels/units and minor textual errors, and refined the presentation and scope of the conclusions derived from the acute EBV infection experiments (including additional clarification of the time window and corresponding immune responses). Below, we provide a detailed point-by-point response to each comment. Reviewers’ comments are reproduced in italics, followed by our responses. All changes are highlighted in the revised manuscript and are indicated by page/line numbers and figure references where applicable.\nReviewers’ Comments to Author\nReviewer #1:\nThis reviewer is overall satisfied with the work that has been done by the authors.\nRegarding the bioinformatics pipeline, the authors have addressed my concerns in great detail, and I am fully satisfied with their responses. The release of the public GitHub repository and the long-term archive on Zenodo provides the transparency and transparency needed for true reproducibility. Regarding the cell type annotation strategy (point b), the two-stage procedure is now clearly explained. Specifically, I find the implementation of a label-probability threshold combined with cluster-level majority voting and subsequent manual curation to be a robust approach that ensures high-confidence cell type assignments. I have no further concerns about the data reproducibility issue, and I believe the manuscript has been significantly improved by these additions.\nResponse: We sincerely thank the reviewer for the positive evaluation of our revised manuscript. We greatly appreciate the reviewer’s constructive suggestions from the previous round, which were very helpful in improving the clarity, transparency, and overall quality of the work. We are pleased that the reviewer finds these additions satisfactory and has no further concerns regarding data reproducibility.\nFor the experimental part, there are still few concerns that need to be addressed:\n(1) Figure R1, the ELISA standard curves for IL-1β and TNF-α range from 0 to 300 pg/mL, whereas the reported results reach values of 15,000 or even 40,000 pg/mL. However, based on the yellow color intensity observed in the plates shown on the left, the measured signals appear to be well below the upper limit of the standard curves. This discrepancy suggests that there may be errors in the data analysis or reporting.\nResponse: Thank you for this important comment. We apologize for the potential confusion. We would like to clarify that all ELISA measurements were performed on appropriately diluted culture supernatants, such that the OD450 values of the diluted samples fell within the linear working range and within the standard-curve range (0–300 pg/mL). The concentrations reported in the manuscript correspond to the calculated values for the original (undiluted) samples, obtained by multiplying the concentration from the standard curve by the dilution factor. Specifically, in Figure R1, samples were diluted 100-fold for IL‑1β, 5-fold for IL‑6 and 200-fold for TNF‑α. In Figure R4B, samples were diluted 5-fold for IL‑1β, and 50-fold for IL‑6 and TNF‑α. To prevent misunderstanding, we have revised the Supplementary Figure legends to explicitly state the dilution factors and clarify that all reported concentrations are dilution-corrected values for the original samples (Supplementary Figure legends, lines 243 and 263).\n(2) In Figure R3, the experimental setup does not specify the duration over which cell proliferation and differentiation were evaluated.\nResponse: Thank you for this comment. We agree that the previous version did not clearly specify the duration over which T‑cell proliferation and effector polarization (“differentiation” in the schematic) were assessed. In the revised manuscript and Supplementary Methods, we have added a detailed experimental timeline. Briefly, purified human CD3⁺ T cells were CFSE‑labelled, stimulated with PHA, and cultured with conditioned supernatants from WT or NR1H3‑OE THP‑1 macrophages. T‑cell proliferation was assessed by CFSE dilution after 72 h (3 days) by flow cytometry. T‑cell effector polarization was assessed by intracellular IFN‑γ staining within CD8⁺ T cells after 24 h of exposure to the conditioned supernatants. In addition, the IFN‑γ ELISPOT assay was performed with a 24‑h incubation period. These clarifications have been incorporated into the updated Supplementary Methods and the figure legend(Supplementary Fig. 17C legends, lines 251).\n(3) CFSE staining typically results in multiple distinct peaks, each corresponding to a successive cell division. However, in the reported figure, the cells appear not to have proliferated in either the treated or untreated conditions. Therefore, the appropriate conclusion would be that no cell proliferation occurred, rather than that no difference in proliferation was observed between the two conditions. The Facs plot in Figure C shows a shift towards the bottom part in the OE treatment. Backgating should be reported to classify those cells as perfectly viable.\nResponse: Thank you for raising this important point. We agree with the reviewer’s assessment. In the CFSE assay, the profiles did not display clear stepwise dilution peaks in either condition, indicating that robust/detectable T‑cell proliferation was not observed under the current assay conditions, rather than supporting a conclusion of “no difference in proliferation” between groups. Accordingly, we have revised the Results text and the corresponding figure legend to reflect the correct interpretation, stating that proliferation was negligible/not detectable in both conditions (Results section Line 273; Supplementary Fig. 17F legend Line 261).\nRegarding the apparent downward shift in the NR1H3‑OE condition, all flow‑cytometry analyses were performed on live singlets using a sequential gating strategy (FSC/SSC → singlets → viability dye → CD3⁺ → CD8⁺ → IFN‑γ). To directly address the reviewer’s request, we provide the full gating strategy together with backgating of the CD8⁺IFN‑γ⁺ events onto the viability dye and FSC/SSC gates in the response letter (Response Fig. R1), confirming that the quantified populations are derived from viable cells. Therefore, the observed shift reflects reduced IFN‑γ staining rather than reduced cell viability.\nRFigure 1\n(4) Rfigure4, have the Elisa data been normalized by cell count?\nResponse: Thank you for raising this important point. For the ELISA experiment in RFigure 4, macrophages were plated at the same seeding density in all wells and treated under identical conditions. At the time of supernatant collection, we assessed cell number and viability, and did not detect a significant difference between the T0901317‑treated and vehicle groups. Therefore, the ELISA values are reported as cytokine concentrations in the culture supernatants (pg/mL) without additional normalization to cell number, as the compared conditions contained comparable numbers of viable cells. To avoid ambiguity, we have clarified this point in the Supplementary Fig. 17H legend (Line 269: “Macrophages were seeded at equal densities; cell number at harvest was comparable between groups; cytokines are shown as pg/mL in conditioned supernatants.”).\n(5) Fig 2 C, log normalized – unit is missing\nResponse: We thank the reviewer for pointing this out. In Fig. 2C, the colour scale reflects the mean log-normalized gene expression values, which are dimensionless and therefore do not have an associated unit. To avoid any confusion, we have clarified this directly in the Fig. 2C legend using the exact wording below: “The colour scale corresponds to the mean log-normalized expression (unitless) of specific genes.”\n(6) Sup F17 cultrue should read culture and “Cultural supernatant” should read “conditioned supernatant”\nResponse: Thank you. We have corrected “Cultural supernatant” with “conditioned supernatant” in Sup F17.\nReviewer #2\nThe manuscript is well revised following my comment and now I agree that the manuscript has reached the level for the publication.\nResponse: We sincerely thank the reviewer for the positive evaluation and for acknowledging that the revised manuscript has reached the level for publication. We greatly appreciate the reviewer’s time and constructive comments, which have helped improve the quality and clarity of our work.\nReviewer #3\nRemarks to the Author:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nRegarding the code availability, the authors have addressed my concerns in great detail, and I am fully satisfied with their responses. The release of the public GitHub repository and the long-term archive on Zenodo provides the transparency and transparency needed for true reproducibility. Regarding the cell type annotation strategy (point b), the two-stage procedure is now clearly explained. Specifically, I find the implementation of a label-probability threshold combined with cluster-level majority voting and subsequent manual curation to be a robust approach that ensures high-confidence cell type assignments. I have no further concerns about the data reproducibility issue, and I believe the manuscript has been significantly improved by these additions.\nResponse: We sincerely thank you for the thoughtful and constructive assessment, as well as for the time invested in co-reviewing this manuscript as part of the Nature Communications initiative. We are pleased that the additions on code availability (GitHub and Zenodo) and the clarified two-stage cell-type annotation strategy have addressed all concerns.\nReviewer #5\nRemarks to the Author:\nIn their revised manuscript the authors have addressed most of this reviewer’s concerns. Namely, that provide more detailed insights into tonsillar architecture, anatomical location and cellular content, which indeed show similarities with human tonsils. Furthermore, they acknowledge that intravenous EBV infection is not physiologic and also compromises investigating viral infection as well as responding immune cells in the tonsils. It seems reasonable that an oropharyngeal route of EBV infection will be investigated in the future. Following this reviewer’s request the authors have increased the number of EBV infected tree shrews following viral infection at 3 and 7 days. Alterations of the gene expression by NR1H3 (LXRa+) macrophages was confirmed. In response to the request for macrophage subset depletion the authors point towards lack of suitable reagents, provide additional in vitro date on NR1H3 manipulation in THP1 cells and now cite published evidence from mouse gamma-herpesvirus infection that is less well controlled without NR1H3. Finally, they provide a T cell analysis, suggesting induction of a starting CTL response by EBV infection of tree shrews. However, the kinetic is with 3 and 7 dpi rather short, considering maximal CD8+ T cell expansion during infectious mononucleosis 4 to 6 weeks after primary infection. In summary, the authors addressed sufficiently major concerns 1, half of 3, 4 and 5, as well as minor concerns 1 and 2 in response to this reviewer. Therefore, the manuscript is significantly improved.\nResponse: We sincerely thank the reviewer for the detailed assessment and for recognizing that the revised manuscript addresses most concerns and has been significantly improved. We appreciate the reviewer’s comments on (i) the expanded description of tonsillar architecture, anatomical location, and cellular composition; (ii) our explicit discussion of the limitations of intravenous EBV inoculation and the rationale for pursuing an oropharyngeal route in future work; and (iii) the increased number of EBV-infected tree shrews analyzed at 3 and 7 dpi, which strengthened our conclusions regarding NR1H3 macrophage-associated transcriptional changes.\nRegarding the request for macrophage subset depletion, we agree that this would be valuable; however, as noted, suitable reagents are currently not available for this model. We therefore complemented the in vivo findings with additional in vitro data on NR1H3 manipulation and by citing relevant in vivo evidence supporting a functional role for NR1H3 in controlling herpesvirus infection.\nFinally, we appreciate the reviewer’s feedback on the T-cell analysis and fully agree that the 3 and 7 dpi window captures only an early phase of infection and is short relative to the peak CD8+ T-cell expansion reported during infectious mononucleosis several weeks after primary infection. To make this limitation explicit, we have revised the Discussion (Lines 423–425) to clarify that our EBV infection model primarily focuses on the acute phase and does not capture longer-term adaptive immune kinetics (including CD8+ T-cell expansion) or NR1H3-mediated regulation during chronic infection. We also note that extended follow-up will be an important direction for future work to resolve longer-term immune dynamics, including later-stage CD8+ T-cell responses.\nChanges in the manuscript: We have added text at Line 423: “Furthermore, our EBV infection model focused primarily on the acute phase, leaving the longer-term kinetics of adaptive immune responses (including CD8+ T-cell expansion) and NR1H3-mediated regulation during chronic infection unexplored.”\nDear Editors,\nWe thank you for the positive editorial decision and we are pleased to submit the revised version of our manuscript. We have carefully addressed all remaining reviewer comments and the journal’s editorial requests. Below we provide a point-by-point response to the reviewers. For clarity, reviewer comments are reproduced in italics, followed by our responses.\nReviewer #1\nComment: “The authors fulfilled all the requests from this reviewer and the manuscript is now ready for publication.”\nResponse: We sincerely thank the you for the careful evaluation of our work and for confirming that all requests have been addressed.\nReviewer #5\nComment: The reviewer notes the value of the tree shrew immune atlas, but highlights that the main weakness of the previous versions concerned the interpretation of the EBV infection data. In the second revision, the reviewer acknowledges that we have clarified the limitation that 3–7 dpi is too early to evaluate CD8⁺ T-cell expansion and that later stages will be addressed in future studies, and considers the manuscript improved.\nResponse: We thank you for the constructive evaluation and for acknowledging the improvements in our revised manuscript. In response to the reviewer’s previous concern regarding the interpretation of the EBV data, we have further clarified in the Discussion that our EBV analysis focuses on the acute 3–7 dpi window and we avoid any claims about later CD8⁺ T-cell expansion kinetics, which will be addressed in future work.\nSincerely,\nWei Xia and Anzhou Tang, on behalf of all authors"}, {"round": 2, "reviewer_comments": {"reviewer_1": "This reviewer is overall satisfied with the work that has been done by the authors.\n\nRegarding the bioinformatics pipeline, the authors have addressed my concerns in great detail, and I am fully satisfied with their responses. The release of the public GitHub repository and the long-term archive on Zenodo provides the transparency and transparency needed for true reproducibility. Regarding the cell type annotation strategy (point b), the two-stage procedure is now clearly explained. Specifically, I find the implementation of a label-probability threshold combined with cluster-level majority voting and subsequent manual curation to be a robust approach that ensures high-confidence cell type assignments. I have no further concerns about the data reproducibility issue, and I believe the manuscript has been significantly improved by these additions.\n\nFor the experimental part, there are still few concerns that need to be addressed:\n\nFigure R1, the ELISA standard curves for IL-1β and TNF-α range from 0 to 300 pg/mL, whereas the reported results reach values of 15,000 or even 40,000 pg/mL. However, based on the yellow color intensity observed in the plates shown on the left, the measured signals appear to be well below the upper limit of the standard curves. This discrepancy suggests that there may be errors in the data analysis or reporting.\n\nIn Figure R3, the experimental setup does not specify the duration over which cell proliferation and differentiation were evaluated.\n\nCFSE staining typically results in multiple distinct peaks, each corresponding to a successive cell division. However, in the reported figure, the cells appear not to have proliferated in either the treated or untreated conditions. Therefore, the appropriate conclusion would be that no cell proliferation occurred, rather than that no difference in proliferation was observed between the two conditions. The Facs plot in Figure C shows a shift towards the bottom part in the OE treatment. Backgating should be reported to classify those cells as perfectly viable.\n\nRfigure4, have the Elisa data been normalized by cell count?\n\nFig 2 C, log normalized – unit is missing\n\nSup F17 cultrue should read culture and “Cultural supernatant” should read “conditioned supernatant”", "reviewer_2": "The manuscript is well revised following my comment and now I agree that the manuscript has reached the level for the publication.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\n\nRegarding the code availability, the authors have addressed my concerns in great detail, and I am fully satisfied with their responses. The release of the public GitHub repository and the long-term archive on Zenodo provides the transparency and transparency needed for true reproducibility. Regarding the cell type annotation strategy (point b), the two-stage procedure is now clearly explained. Specifically, I find the implementation of a label-probability threshold combined with cluster-level majority voting and subsequent manual curation to be a robust approach that ensures high-confidence cell type assignments. I have no further concerns about the data reproducibility issue, and I believe the manuscript has been significantly improved by these additions.", "reviewer_5": "This new study describes a single cell immune atlas for the tree shrew (Tupaia belangeri). Tree shrews are an emerging small animal model that are evolutionarily closer to primates than mice. As such, tree shrews represent a useful, viable alternative model for infectious disease and aging studies that may be otherwise challenging to carry out in murine or non-human primate model systems. The investigators profiled single cells from 12 different tissues of 3 adult tree shrews, identifying 39 different cell types. Notably, bulk RNAseq studies have been performed in spleen, PBMCs, and other tissues from the tree shrew, and single cell transcriptome atlases have been recently published for brain and neurological tissues. A strength of this study is that it represents the first tree shrew immune-focused single cell atlas which can provide an important platform for understanding systemic transcriptional responses to infectious diseases in future studies. Aside from the substantial efforts that went into generating this resource, there are some key weaknesses and major concerns with the overall conclusions derived from the interpretation of the EBV infection data.\n\nIn their revised manuscript the authors have addressed most of this reviewer’s concerns. Namely, that provide more detailed insights into tonsillar architecture, anatomical location and cellular content, which indeed show similarities with human tonsils. Furthermore, they acknowledge that intravenous EBV infection is not physiologic and also compromises investigating viral infection as well as responding immune cells in the tonsils. It seems reasonable that an oropharyngeal route of EBV infection will be investigated in the future. Following this reviewer’s request the authors have increased the number of EBV infected tree shrews following viral infection at 3 and 7 days. Alterations of the gene expression by NR1H3 (LXRa+) macrophages was confirmed. In response to the request for macrophage subset depletion the authors point towards lack of suitable reagents, provide additional in vitro date on NR1H3 manipulation in THP1 cells and now cite published evidence from mouse gamma-herpesvirus infection that is less well controlled without NR1H3. Finally, they provide a T cell analysis, suggesting induction of a starting CTL response by EBV infection of tree shrews. However, the kinetic is with 3 and 7 dpi rather short, considering maximal CD8+ T cell expansion during infectious mononucleosis 4 to 6 weeks after primary infection. In summary, the authors addressed sufficiently major concerns 1, half of 3, 4 and 5, as well as minor concerns 1 and 2 in response to this reviewer. Therefore, the manuscript is significantly improved."}, "author_response": "Response to Reviewers’ Comments\nWe thank the Editor and the reviewers for their careful evaluation of our manuscript and for the constructive comments provided in this round. We have revised the manuscript accordingly and addressed all remaining points raised by the reviewers. In particular, we have clarified the experimental design and data interpretation for the ELISA and CFSE assays, corrected figure labels/units and minor textual errors, and refined the presentation and scope of the conclusions derived from the acute EBV infection experiments (including additional clarification of the time window and corresponding immune responses). Below, we provide a detailed point-by-point response to each comment. Reviewers’ comments are reproduced in italics, followed by our responses. All changes are highlighted in the revised manuscript and are indicated by page/line numbers and figure references where applicable.\n\nReviewers’ Comments to Author\nReviewer #1:\nThis reviewer is overall satisfied with the work that has been done by the authors.\nRegarding the bioinformatics pipeline, the authors have addressed my concerns in great detail, and I am fully satisfied with their responses. The release of the public GitHub repository and the long-term archive on Zenodo provides the transparency and transparency needed for true reproducibility. Regarding the cell type annotation strategy (point b), the two-stage procedure is now clearly explained. Specifically, I find the implementation of a label-probability threshold combined with cluster-level majority voting and subsequent manual curation to be a robust approach that ensures high-confidence cell type assignments. I have no further concerns about the data reproducibility issue, and I believe the manuscript has been significantly improved by these additions.\nResponse: We sincerely thank the reviewer for the positive evaluation of our revised manuscript. We greatly appreciate the reviewer’s constructive suggestions from the previous round, which were very helpful in improving the clarity, transparency, and overall quality of the work. We are pleased that the reviewer finds these additions satisfactory and has no further concerns regarding data reproducibility.\nFor the experimental part, there are still few concerns that need to be addressed:\n(1) Figure R1, the ELISA standard curves for IL-1β and TNF-α range from 0 to 300 pg/mL, whereas the reported results reach values of 15,000 or even 40,000 pg/mL. However, based on the yellow color intensity observed in the plates shown on the left, the measured signals appear to be well below the upper limit of the standard curves. This discrepancy suggests that there may be errors in the data analysis or reporting.\nResponse: Thank you for this important comment. We apologize for the potential confusion. We would like to clarify that all ELISA measurements were performed on appropriately diluted culture supernatants, such that the OD450 values of the diluted samples fell within the linear working range and within the standard-curve range (0–300 pg/mL). The concentrations reported in the manuscript correspond to the calculated values for the original (undiluted) samples, obtained by multiplying the concentration from the standard curve by the dilution factor. Specifically, in Figure R1, samples were diluted 100-fold for IL‑1β, 5-fold for IL‑6 and 200-fold for TNF‑α. In Figure R4B, samples were diluted 5-fold for IL‑1β, and 50-fold for IL‑6 and TNF‑α. To prevent misunderstanding, we have revised the Supplementary Figure legends to explicitly state the dilution factors and clarify that all reported concentrations are dilution-corrected values for the original samples (Supplementary Figure legends, lines 243 and 263).\n(2) In Figure R3, the experimental setup does not specify the duration over which cell proliferation and differentiation were evaluated.\nResponse: Thank you for this comment. We agree that the previous version did not clearly specify the duration over which T‑cell proliferation and effector polarization (“differentiation” in the schematic) were assessed. In the revised manuscript and Supplementary Methods, we have added a detailed experimental timeline. Briefly, purified human CD3⁺ T cells were CFSE‑labelled, stimulated with PHA, and cultured with conditioned supernatants from WT or NR1H3‑OE THP‑1 macrophages. T‑cell proliferation was assessed by CFSE dilution after 72 h (3 days) by flow cytometry. T‑cell effector polarization was assessed by intracellular IFN‑γ staining within CD8⁺ T cells after 24 h of exposure to the conditioned supernatants. In addition, the IFN‑γ ELISPOT assay was performed with a 24‑h incubation period. These clarifications have been incorporated into the updated Supplementary Methods and the figure legend(Supplementary Fig. 17C legends, lines 251).\n(3) CFSE staining typically results in multiple distinct peaks, each corresponding to a successive cell division. However, in the reported figure, the cells appear not to have proliferated in either the treated or untreated conditions. Therefore, the appropriate conclusion would be that no cell proliferation occurred, rather than that no difference in proliferation was observed between the two conditions. The Facs plot in Figure C shows a shift towards the bottom part in the OE treatment. Backgating should be reported to classify those cells as perfectly viable.\nResponse: Thank you for raising this important point. We agree with the reviewer’s assessment. In the CFSE assay, the profiles did not display clear stepwise dilution peaks in either condition, indicating that robust/detectable T‑cell proliferation was not observed under the current assay conditions, rather than supporting a conclusion of “no difference in proliferation” between groups. Accordingly, we have revised the Results text and the corresponding figure legend to reflect the correct interpretation, stating that proliferation was negligible/not detectable in both conditions (Results section Line 273; Supplementary Fig. 17F legend Line 261).\nRegarding the apparent downward shift in the NR1H3‑OE condition, all flow‑cytometry analyses were performed on live singlets using a sequential gating strategy (FSC/SSC → singlets → viability dye → CD3⁺ → CD8⁺ → IFN‑γ). To directly address the reviewer’s request, we provide the full gating strategy together with backgating of the CD8⁺IFN‑γ⁺ events onto the viability dye and FSC/SSC gates in the response letter (Response Fig. R1), confirming that the quantified populations are derived from viable cells. Therefore, the observed shift reflects reduced IFN‑γ staining rather than reduced cell viability.\nRFigure 1\n(4) Rfigure4, have the Elisa data been normalized by cell count?\nResponse: Thank you for raising this important point. For the ELISA experiment in RFigure 4, macrophages were plated at the same seeding density in all wells and treated under identical conditions. At the time of supernatant collection, we assessed cell number and viability, and did not detect a significant difference between the T0901317‑treated and vehicle groups. Therefore, the ELISA values are reported as cytokine concentrations in the culture supernatants (pg/mL) without additional normalization to cell number, as the compared conditions contained comparable numbers of viable cells. To avoid ambiguity, we have clarified this point in the Supplementary Fig. 17H legend (Line 269: “Macrophages were seeded at equal densities; cell number at harvest was comparable between groups; cytokines are shown as pg/mL in conditioned supernatants.”).\n(5) Fig 2 C, log normalized – unit is missing\nResponse: We thank the reviewer for pointing this out. In Fig. 2C, the colour scale reflects the mean log-normalized gene expression values, which are dimensionless and therefore do not have an associated unit. To avoid any confusion, we have clarified this directly in the Fig. 2C legend using the exact wording below: “The colour scale corresponds to the mean log-normalized expression (unitless) of specific genes.”\n(6) Sup F17 cultrue should read culture and “Cultural supernatant” should read “conditioned supernatant”\nResponse: Thank you. We have corrected “Cultural supernatant” with “conditioned supernatant” in Sup F17.\nReviewer #2\nThe manuscript is well revised following my comment and now I agree that the manuscript has reached the level for the publication.\nResponse: We sincerely thank the reviewer for the positive evaluation and for acknowledging that the revised manuscript has reached the level for publication. We greatly appreciate the reviewer’s time and constructive comments, which have helped improve the quality and clarity of our work.\nReviewer #3\nRemarks to the Author:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nRegarding the code availability, the authors have addressed my concerns in great detail, and I am fully satisfied with their responses. The release of the public GitHub repository and the long-term archive on Zenodo provides the transparency and transparency needed for true reproducibility. Regarding the cell type annotation strategy (point b), the two-stage procedure is now clearly explained. Specifically, I find the implementation of a label-probability threshold combined with cluster-level majority voting and subsequent manual curation to be a robust approach that ensures high-confidence cell type assignments. I have no further concerns about the data reproducibility issue, and I believe the manuscript has been significantly improved by these additions.\nResponse: We sincerely thank you for the thoughtful and constructive assessment, as well as for the time invested in co-reviewing this manuscript as part of the Nature Communications initiative. We are pleased that the additions on code availability (GitHub and Zenodo) and the clarified two-stage cell-type annotation strategy have addressed all concerns.\nReviewer #5\nRemarks to the Author:\nIn their revised manuscript the authors have addressed most of this reviewer’s concerns. Namely, that provide more detailed insights into tonsillar architecture, anatomical location and cellular content, which indeed show similarities with human tonsils. Furthermore, they acknowledge that intravenous EBV infection is not physiologic and also compromises investigating viral infection as well as responding immune cells in the tonsils. It seems reasonable that an oropharyngeal route of EBV infection will be investigated in the future. Following this reviewer’s request the authors have increased the number of EBV infected tree shrews following viral infection at 3 and 7 days. Alterations of the gene expression by NR1H3 (LXRa+) macrophages was confirmed. In response to the request for macrophage subset depletion the authors point towards lack of suitable reagents, provide additional in vitro date on NR1H3 manipulation in THP1 cells and now cite published evidence from mouse gamma-herpesvirus infection that is less well controlled without NR1H3. Finally, they provide a T cell analysis, suggesting induction of a starting CTL response by EBV infection of tree shrews. However, the kinetic is with 3 and 7 dpi rather short, considering maximal CD8+ T cell expansion during infectious mononucleosis 4 to 6 weeks after primary infection. In summary, the authors addressed sufficiently major concerns 1, half of 3, 4 and 5, as well as minor concerns 1 and 2 in response to this reviewer. Therefore, the manuscript is significantly improved.\nResponse: We sincerely thank the reviewer for the detailed assessment and for recognizing that the revised manuscript addresses most concerns and has been significantly improved. We appreciate the reviewer’s comments on (i) the expanded description of tonsillar architecture, anatomical location, and cellular composition; (ii) our explicit discussion of the limitations of intravenous EBV inoculation and the rationale for pursuing an oropharyngeal route in future work; and (iii) the increased number of EBV-infected tree shrews analyzed at 3 and 7 dpi, which strengthened our conclusions regarding NR1H3 macrophage-associated transcriptional changes.\nRegarding the request for macrophage subset depletion, we agree that this would be valuable; however, as noted, suitable reagents are currently not available for this model. We therefore complemented the in vivo findings with additional in vitro data on NR1H3 manipulation and by citing relevant in vivo evidence supporting a functional role for NR1H3 in controlling herpesvirus infection.\nFinally, we appreciate the reviewer’s feedback on the T-cell analysis and fully agree that the 3 and 7 dpi window captures only an early phase of infection and is short relative to the peak CD8+ T-cell expansion reported during infectious mononucleosis several weeks after primary infection. To make this limitation explicit, we have revised the Discussion (Lines 423–425) to clarify that our EBV infection model primarily focuses on the acute phase and does not capture longer-term adaptive immune kinetics (including CD8+ T-cell expansion) or NR1H3-mediated regulation during chronic infection. We also note that extended follow-up will be an important direction for future work to resolve longer-term immune dynamics, including later-stage CD8+ T-cell responses.\nChanges in the manuscript: We have added text at Line 423: “Furthermore, our EBV infection model focused primarily on the acute phase, leaving the longer-term kinetics of adaptive immune responses (including CD8+ T-cell expansion) and NR1H3-mediated regulation during chronic infection unexplored.”\nDear Editors,\nWe thank you for the positive editorial decision and we are pleased to submit the revised version of our manuscript. We have carefully addressed all remaining reviewer comments and the journal’s editorial requests. Below we provide a point-by-point response to the reviewers. For clarity, reviewer comments are reproduced in italics, followed by our responses.\nReviewer #1\nComment: “The authors fulfilled all the requests from this reviewer and the manuscript is now ready for publication.”\nResponse: We sincerely thank the you for the careful evaluation of our work and for confirming that all requests have been addressed.\nReviewer #5\nComment: The reviewer notes the value of the tree shrew immune atlas, but highlights that the main weakness of the previous versions concerned the interpretation of the EBV infection data. In the second revision, the reviewer acknowledges that we have clarified the limitation that 3–7 dpi is too early to evaluate CD8⁺ T-cell expansion and that later stages will be addressed in future studies, and considers the manuscript improved.\nResponse: We thank you for the constructive evaluation and for acknowledging the improvements in our revised manuscript. In response to the reviewer’s previous concern regarding the interpretation of the EBV data, we have further clarified in the Discussion that our EBV analysis focuses on the acute 3–7 dpi window and we avoid any claims about later CD8⁺ T-cell expansion kinetics, which will be addressed in future work.\nSincerely,\nWei Xia and Anzhou Tang, on behalf of all authors"}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors fulfilled all the requests from this reviewer and the manuscript is now ready for publication.", "reviewer_5": "The reviewer notes the value of the tree shrew immune atlas, but highlights that the main weakness of the previous versions concerned the interpretation of the EBV infection data. In the second revision, the reviewer acknowledges that we have clarified the limitation that 3–7 dpi is too early to evaluate CD8⁺ T-cell expansion and that later stages will be addressed in future studies, and considers the manuscript improved."}, "author_response": "Dear Editors,\nWe thank you for the positive editorial decision and we are pleased to submit the revised version of our manuscript. We have carefully addressed all remaining reviewer comments and the journal’s editorial requests. Below we provide a point-by-point response to the reviewers. For clarity, reviewer comments are reproduced in italics, followed by our responses.\nReviewer #1\nComment: “The authors fulfilled all the requests from this reviewer and the manuscript is now ready for publication.”\nResponse: We sincerely thank the you for the careful evaluation of our work and for confirming that all requests have been addressed.\nReviewer #5\nComment: The reviewer notes the value of the tree shrew immune atlas, but highlights that the main weakness of the previous versions concerned the interpretation of the EBV infection data. In the second revision, the reviewer acknowledges that we have clarified the limitation that 3–7 dpi is too early to evaluate CD8⁺ T-cell expansion and that later stages will be addressed in future studies, and considers the manuscript improved.\nResponse: We thank you for the constructive evaluation and for acknowledging the improvements in our revised manuscript. In response to the reviewer’s previous concern regarding the interpretation of the EBV data, we have further clarified in the Discussion that our EBV analysis focuses on the acute 3–7 dpi window and we avoid any claims about later CD8⁺ T-cell expansion kinetics, which will be addressed in future work.\nSincerely,\nWei Xia and Anzhou Tang, on behalf of all authors"}]

Decision letter