Animal husbandry
All four grey mouse lemurs (M. murinus) included in this atlas originated from the closed captive breeding colony at the Muséum National d’Histoire Naturelle in Brunoy, France, and were transferred together to the University of Texas (Austin) in 2009 and then to Stanford University in 2015 and maintained for non-invasive phenotyping and genetic research as approved by the Stanford University Administrative Panel on Laboratory Animal Care (APLAC number 27439) and in accordance with the Guide for the Care and Use of Laboratory Animals, as previously detailed15. In brief, mouse lemurs were housed indoors in modified marmoset cages with multiple PVC perches and nest boxes in a facility credited by The Association for Assessment and Accreditation of Laboratory Animal Care in a temperature (23.3–24.4 °C) and light-controlled environment (daily 14:10 h and 10:14 h light–dark, alternating every 6 months to synchronize seasonal breeding behaviour and metabolic changes) and were fed ad libitum with fresh fruits and vegetables, crushed primate chow (Teklad Global 20% Protein Primate Diet, 2050, Envigo) and live insect larvae as enrichment items. Animals were socially housed in single-sex groups or individually housed owing to behavioural incompatibility or health management requirements. Health and welfare were monitored daily and clinical care was provided by the Veterinary Service Center, Stanford University, including diagnosis and treatment of spontaneously occurring health conditions. Animals in declining health despite medical care were euthanized for humane reasons as determined by a veterinarian. Before euthanasia, it happened that all four lemurs were living in summer-like long days (14:10 h) for at least 3 months (range 3–6 months), and all showed standard activity patterns without signs of torpor. Given these individuals were housed at constant temperature conditions and fed a non-calorie-restrictive diet, spontaneous torpor was not observed in any of the analysed lemurs throughout their time in the Stanford colony.
Tissue collection and processing
Animals in declining health who did not respond to standard therapy were euthanized by pentobarbital overdose under isoflurane anaesthesia as previously described15. Before euthanasia, a veterinary examination was performed, and animal body weight and electrocardiogram data were obtained (KardiaMobile 6L, AliveCor). Blood was immediately collected by cardiocentesis for serum chemistry, complete blood count, biobanking and scRNA-seq. In three animals (L2–L4), transcardial perfusion of the lungs with PBS was done to reduce circulating cells. Organs and tissues were sequentially removed and divided by a veterinary pathologist. One sample of each tissue was immediately placed in formalin fixative for histopathology15, and a second was embedded in optimal cutting temperature compound and then flash-frozen on dry ice and stored at −80 °C for biobanking. A third sample was placed directly in cold (4 °C) PBS pH 7.4 and immediately distributed to the tissue expert for cell dissociation and preparation for scRNA-seq as detailed below. Additional diagnostics, such as microbiological cultures, were performed where clinically indicated. The entire necropsy was completed within 1–2 h, with ischaemia-sensitive tissues prioritized as described in the Supplementary Methods.
Histological and pathological analysis
Tissues were immersion-fixed in 10% neutral-buffered formalin for 72 h. Formalin-fixed tissues were processed routinely, embedded in paraffin, sectioned at 5 µm and stained with haematoxylin and eosin (H&E). The following tissues were analysed: heart, aorta, lungs, trachea, thyroid gland, parathyroid gland, kidneys, urinary bladder, male reproductive tract (testicle, epididymis, seminal vesicle, prostate and penile urethra), female reproductive tract (uterus, cervix, vagina and ovaries), salivary glands, tongue, epiglottis, oesophagus, stomach, small and large intestine, liver (with gallbladder), adrenal gland, spleen, lymph nodes, white adipose, brown adipose, bone, spinal cord, eyes and bone marrow. Selected tissues were stained with Von Kossa (for mineralization), Masson’s trichrome (for collagen), Congo Red (for amyloid) and Gram stain (for bacteria) as part of the pathological analysis. H&E-stained slides were scanned with a Leica Aperio AT2 high-volume digital whole-slide scanner (×40 objective), uploaded into Napari image viewer44, software adapted by CZB and posted on the Tabula Microcebus portal.
Preparation of single-cell suspensions and FACS for scRNA-seq
Fresh tissue samples obtained as described above were placed on ice, delivered to tissue experts and immediately dissociated and processed into single-cell suspensions, except for samples from L3, which were kept cool overnight after necropsy and processed the next morning (Supplementary Methods). For each solid tissue, this process involved a standard combination of enzymatic digestion and mechanical disruption methods that were optimized for the specific tissue, many of which were adapted from procedures used for the corresponding mouse tissue12,13. For blood, immune cells were isolated using a high-density Ficoll gradient (Histoplaque-1119, Sigma-Aldrich) to include peripheral blood mononuclear cells and polymorphonuclear leukocytes14.
The specific protocols for each of the 27 tissues are detailed in the Supplementary Methods. The cell number and concentration for each single-cell suspension were determined by manual counting using a haemocytometer and then adjusted with 2% FBS in PBS to a target concentration of about 106 cells per ml. Samples were then used for droplet-based 10x library preparation and/or flow sorted for single live cells (Sytox blue negative; ThermoFisher, S34857) for plate-based SS2 library preparation (Supplementary Fig. 5). To enrich for cardiomyocytes, the standard procedure for cardiac cell isolation was supplemented by hand-picking cardiomyocytes (Supplementary Methods). Residual cell suspensions were diluted 1:1 with serum-free Bambanker cell freezing medium (GC Lymphotec, BB01) and cryopreserved at −80 °C.
scRNA-seq library preparation, quality control and sequencing
For 10x, single cells were profiled using the 10x Genomics scRNA-seq pipeline (Chromium Single Cell 3′ Library and Gel Bead v.2 Chemistry kit) and sequenced on a NovaSeq 6000 System as previously described12,13,14 and detailed in the Supplementary Methods. For SS2, single cells were sorted into 384-well or 96-well lysis plates, reverse transcribed to cDNA and amplified, as previously described12,13. cDNA libraries were prepared using a Nextera XT Library Sample Preparation kit (Illumina, FC-131-1096) or (for L4) an in-house protocol detailed in the Supplementary Methods; no significant differences between protocols were observed in library read depth or quality. Pooling of individual libraries and subsequent quality control and DNA sequencing were done as previously described12,13,14 with minor modifications (Supplementary Methods). Both 10x and SS2 libraries were sequenced to achieve saturation on an Illumina NovaSeq 6000 system (10x, 26 bp and 90 bp paired-end reads; SS2, 2 × 100 bp paired-end reads).
Genome alignment of scRNA-seq reads and gene counts
The M. murinus genome assembly (Mmur 3.0, accession: GCF_000165445.2; annotation: NCBI Refseq annotation release 101) with NCBI annotation release 101 (date acquired, 21 September 2018) was used for downstream alignment and data analyses. A total of 31,509 genes were detected, including annotated genes and unannotated loci but excluding mitochondrial and Y chromosome genes (unannotated at our acquisition date).
For 10x samples, downstream data were processed using standard methods with Cell Ranger (v.2.2, 10x Genomics). Raw base call files directly generated by the NovaSeq instrument were demultiplexed and converted to fastq files and then aligned to the 10x genome index, with barcode and UMI counting performed to generate a gene counts table. Alignment files were outputted in standard BAM format.
For SS2 samples, demultiplexed fastq files were mapped to the genome using STAR aligner (v.2.6.1a). In brief, the genome FASTA file was augmented with ERCC sequences to create a STAR genome index with 99 bp overhangs (optimized for Illumina 2 × 100 bp paired-end reads). Two-pass mapping was executed, in which the first pass identified splice junctions that were added to the gene reference to improve second pass mapping, with specific STAR options and parameters detailed in the Supplementary Methods.
Contamination filtering of 10x data
We performed stringent contamination filtering to resolve cross-sample contamination in an Illumina sequencing run caused by cell barcode hopping among multiplexed 10x samples45,46. Such cross-sample contamination can occur when low levels of ambient mRNA containing the 10x cell barcode in one sample gets added onto the transcript of other samples during Illumina sequencing amplification, which results in the incorrect assignment of a cell barcode to other samples. Hence, in subsequent analyses, a cell from one tissue could falsely appear as multiple cells from different tissues (or samples). To exclude such artefacts, for each sequencing run, we identified all cell barcodes that were assigned to multiple samples. Then for each such barcode identified, we compared the number of UMIs in each sample. If there was one dominant sample index (that is, the number of UMIs of the dominant sample was ten times or more greater than that of the second most abundant sample), then the cell with the dominant sample index was kept (but labelled in its metadata as ‘potentially contaminated’), whereas all other instances of that ‘cell’ were removed. If there was no dominant sample index, then all instances of the ‘cell’ with that barcode were removed from the dataset. Contamination was not an issue for SS2 samples because they were sequenced using dual unique indices for each cell.
Cell clustering, annotation and cluster markers from scRNA-seq profiles
Cell clustering and annotation of each tissue processed by 10x
Transcriptomic profiles of cells from each tissue and from each individual lemur were clustered separately using Seurat software (v.2.3.0) for R studio (v.3.6.1). We included in this step all cells with >100 genes or >1,000 UMIs detected, a minimal threshold that was used to ensure the inclusion of all cell types, including ones in which the cells (or RNA) were unstable (see below for more stringent criteria used for final cell quality control). For each cell, expression of a gene g was normalized in 10x data as follows: ln(UMIg/UMItotal × 1 × 104 + 1), abbreviated as ln(UP10K + 1); in SS2: ln(readsg/readstotal × 1 × 104 + 1), abbreviated as ln(CP10K + 1). Next, data scaling, dimensionality reduction (PCA), clustering and visualization (t-SNE and UMAP) were performed following the standard Seurat pipeline as previously described14, with parameters including the numbers of principal components (PCs), perplexity and resolution manually adjusted for each iteration of cell clustering. Resultant cell clusters were manually assigned to a compartment (endothelial, epithelial, stromal, lymphoid, myeloid, megakaryocyte–erythroid, haematopoietic (for precursors), neural and germ) on the basis of the expression of the mouse lemur orthologues of canonical marker genes for each compartment in human and mouse (Supplementary Table 1). Clusters that expressed markers from more than one compartment were annotated as ‘doublets’. Cells in each assigned compartment were then separately subclustered, repeating the data processing steps above. To annotate (determine) the cell type of each cluster in a compartment, a list of canonical human and mouse gene markers for each cell type in each tissue was curated from the literature (Supplementary Table 1), including genes previously validated by in situ hybridization and/or immunohistochemistry as well as DEGs selected from recent scRNA-seq studies. The orthologous mouse lemur genes were identified and their expression visualized on the t-SNE plots. On the basis of the enriched expression of marker genes, each cluster of cells in a compartment was manually assigned a cell-type identity. Clusters that contained more than one cell type were further subclustered to better resolve the cell types. Cell types represented by only a small number of cells that did not form a separate cluster were manually curated, aided by the cellxgene gene expression visualization tool47 as detailed below.
Each cluster was assigned both a ‘cell ontology’ cell-type designation (name) using the standardized and structured nomenclature19 and a ‘free annotation’ that resolved biologically significant clusters not contained in the current cell ontology. Free annotations were assigned as follows. In cases when a smaller cluster stemmed off a larger (main) cluster in the t-SNE-embedded space, the smaller cluster was distinguished with one or more DEGs added to the cell-type name (for example, B cell (SOX5+) clustered near the main population of B cells in the pancreas). DEGs driving the subtype clustering were ascertained by Wilcoxon rank-sum tests. In cases when two approximately equal-sized clusters separated on the t-SNE plot, a marker gene was added to the cell-type name for both clusters. Clusters with a small number of cells that contained more than one cell type but could not be partitioned into separate clusters by subclustering with the Louvain algorithm or manually with cellxgene (see the section ‘Integration of datasets across individuals’) were labelled as a ‘mixed’ cell type (for example, the cluster labelled ‘endothelial cell’ in the uterus contains a mixture of artery, vein and capillary cells). Clusters with cells that expressed markers for more than one cell type and it was biologically plausible that they were not a technical artefact (for example, a doublet of two distinct cell types) were labelled as a ‘hybrid’ cell type (for example, the cluster labelled as ‘monocyte–macrophage’ in the trachea contains cells that expressed markers of both cell types and could not be further distinguished based on current molecular definitions of these cell types). After examining the human and mouse markers for all known cell types in a tissue, clusters that could not be assigned a cell type were labelled ‘unknown’, with the tissue, compartment and one or more DEGs added to the cell type name (for example, ‘unknown bone stromal G1 (NGFR+TNNT2+)’ are bone stromal cells that do not correspond to any extant stromal cell type reported for humans or mice). To detect the DEGs of an unknown cell type, we compared the unknown cell type to all other cells of the same compartment and tissue (Extended Data Fig. 4 and Supplementary Fig. 3). Clusters containing a majority of cells that expressed cell proliferation markers (for example, TOP2A, MKI67 and STMN1) were appended the abbreviation ‘PF’. Clusters that separated from a main cluster but did not express any distinguishing markers (other than tRNAs, rRNAs and/or immediate-early genes) and differed only in parameters of technical quality (that is, fewer genes and counts detected per cell) were considered low quality and ‘LQ’ was appended to the cell-type name.
After annotations were assigned, the cut-off for the minimum number of genes per cell was increased from 100 to 500, and only the qualifying cells were further analysed. For most tissues, this more stringent cut-off value only resulted in removal of some erythrocytes and neutrophils. The only exception were cardiac cardiomyocytes, most of which expressed fewer than 500 genes per cell; therefore, separate filtering criteria were applied (Supplementary Methods).
Annotation of each tissue processed by SS2
Cells processed using the SS2 protocol with <500 genes or <5,000 reads were excluded from further analysis, and gene expression levels in the remaining cells were scaled and log transformed as described above for the 10x datasets. Cells from a particular tissue and individual were integrated with the 10x dataset of the same tissue and individual into the same UMAP-embedded space using the FIRM algorithm (detailed below). Cells from SS2 were automatically annotated with the same label as the nearest neighbouring 10x cell. Annotations were manually verified as described in the section ‘Integration of datasets across individuals’, aided by cellxgene gene expression visualization. SS2 datasets for which there were no corresponding 10x dataset from the same individual or tissue were manually annotated using the method described in the section ‘Cell clustering and annotation of each tissue processed by 10x’ for 10x datasets.
Integration of datasets across individuals
For each tissue, the combined 10x and SS2 datasets from each individual were further integrated into the same UMAP-embedded space using the FIRM algorithm18. This step resulted in 27 separate tissue UMAPs, each containing data from up to 4 individuals. To ensure consistency of cell-type labelling across all individuals, annotations were verified or manually adjusted using cellxgene, an interactive tool to visualize and annotate scRNA-seq data (https://chanzuckerberg.github.io/cellxgene/)47.
Integration of datasets across tissues
All 27 tissue-level objects were integrated into a single UMAP-embedded space using the FIRM algorithm. As described above, annotations were verified or manually adjusted in cellxgene to ensure consistency of cell designations across all tissues. In most instances, cells from the same cell type clustered together, irrespective of the tissue of origin, and the same designation was used across all tissues. Occasionally, similar cells types (for example, fibroblasts and macrophages) clustered separately by tissue of origin, which made it challenging to distinguish whether the separation was due to a tissue-level batch effect or because of true biological differences. In these cases, the original tissue-level annotation label was kept for each cluster. In total, 256 cell designations were assigned across the integrated atlas, which, when distinguished by organ of origin (for example, lung versus bladder artery cells), resulted in a total of 768 molecular cell types.
Detection of DEGs for each cell type
We calculated the top 300 DEGs (adjusted P < 0.05) for each cell type in the 10x dataset (represented by at least 5 individual cells after removing doublets, low-quality cells and mixed cell types) using two-tailed Wilcoxon rank-sum tests with Benjamini–Hochberg false discovery rate correction (Supplementary Table 3). We compared each cell type to the following: (1) all other cell types from the same tissue (for example, lung capillary cells compared with all other lung cells; ‘tissue-wide’ comparison); (2) all other cell types from the same compartment of that tissue (lung capillary cells compared with all other lung endothelial cells; ‘tissue-compartment-wide’ comparison); (3) all other cell types from the atlas (lung capillary cells with all other cells in the atlas; ‘atlas-wide’ comparison); and (4) all other cell types from the same compartment across the atlas (lung capillary cells with all other endothelial cells in the atlas; ‘atlas-compartment-wide’ comparison).
FIRM integration
FIRM is a newly developed algorithm that integrates multiple scRNA-seq datasets18 (for example, from different sequencing platforms, tissue types and experimental batches). In brief, FIRM optimizes dataset integration by harmonizing differences in cell-type composition and computing the dataset-specific scale factors for gene-level normalization. Different datasets generally have varied cell-type compositions, which results in dataset disparity when scaling the gene expression levels to the unit variance for each dataset. Different from classical scaling procedures, FIRM computes the scale factors based on subsets of cells that have matched cell-type compositions between datasets. To construct these subsets, FIRM detects paired clusters between datasets based on similar overall gene expression levels and then samples the cells so that paired cell types have the same proportional representation in each dataset. The parameters used for integration are given in the Supplementary Methods. The integrated datasets generated using FIRM showed accurate mixing of shared cell-type identities and preserved the structure of the original datasets, as confirmed by expert manual inspection during cell annotation.
Trajectory analysis
We used two independent methods to characterize spatial and developmental pseudotime cell trajectories: a custom in-house program in Matlab and Slingshot48. For the mouse lemur kidney nephron spatial trajectory, all kidney epithelial cells were included in the analysis, with the exception of podocytes, macula densa cells, intercalated cells and urothelial cells which clustered separately. For the vasa recta endothelium spatial trajectory, all four vasa recta cell types were used. For the spermatogenesis pseudotime trajectory, all seven sperm and sperm progenitor cell types were used. For the myeloid cell developmental pseudotime trajectory, haematopoietic precursor cells and all myeloid cell types except DCs (which did not form part of the continuum) were used. Analysis was performed independently for each trajectory using values from the 10x scRNA-seq profiles of the indicated cells (low-quality cells and technical doublets were excluded) that had been pre-processed (scaled: ln(UP10K + 1), normalized) as described above.
PCA with highly variable genes (dispersion > 0.5) was done with the PCA function of Matlab, and the high-quality PCs (not driven by extreme outlier data points or immediate-early genes) were selected from the top 20 PCs and used to generate a 2D UMAP using cell–cell Euclidean distances as input (https://www.mathworks.com/matlabcentral/fileexchange/71902). The trajectory of the cell continuum was detected as the probability density ridge of the data points in the UMAP, using automated image processing (Matlab Image Processing Toolbox). Any interruptions in the detected density ridge line were manually connected along the direction of the ridge line and guided by previous knowledge of the biological process. The direction of the trajectory was assigned on the basis of expression of marker genes. Individual cells were then aligned to the trajectory by the shortest connecting point to the trajectory; if the trajectory branched (for example, in myeloid cell development), cells were assigned to the closest branch. Individual cells that were too distant from the trajectory (adaptive thresholding along the trajectory) were deemed outliers and removed from further analysis.
To detect genes for which expression followed the trajectory, we calculated Spearman correlation coefficients and corresponding P values (Bonferroni corrected) between the expression level of each gene and 20 preassigned unimodal patterns that smoothly changed along the trajectory (with their single peaks uniformly distributed from the beginning of the trajectory to its end point). Expression patterns of the top ranking (top 1,000 with P < 0.01) and highly variable (dispersion > 0.5) genes were smoothed with a moving average filter and clustered by k-means clustering to detect the major trajectory-dependent expression patterns. The trajectory DEGs were then ranked by the associated cluster (ranked by trajectory location of peak expression), in the cluster by P value from the smallest to the largest, and with the same P value by mean expression level from the highest to the lowest.
For the myeloid cell analysis, four trajectories were independently detected: (1) from haematopoietic precursors to granulocyte–monocyte progenitors; (2) from granulocyte–monocyte progenitors to proliferating and then maturing neutrophils; (3) from granulocyte–monocyte progenitors to proliferating and then maturing monocytes and macrophages; and (4) from megakaryocyte and erythroid progenitors to proliferating and maturing erythroid lineage cells. On the UMAP, trajectory 1 branched into trajectories 2 and 3, so two longer trajectories were generated (1 + 2 and 1 + 3). Differential gene expression analysis was then independently performed for each of the constituent trajectories (1 + 2, neutrophil lineage; 1 + 3, monocyte–macrophage lineage; 4, erythrocyte lineage).
As an alternative method, we applied the Slingshot method48, which first computes the global lineage structure by constructing a cluster-based minimum spanning tree followed by pseudotime inference using simultaneous principal curves to fit smooth branching curves to these lineages. We used the annotated clusters and UMAP coordinates to first obtain a global lineage structure with getLineages, then constructed smooth curves, ordered cells along the trajectory and generated pseudotime values using getCurves. For each tissue, the longest trajectory that incorporated the most clusters was used. For the immune cell trajectories, the neutrophil cluster was subdivided into higher resolution clusters that were then combined to facilitate building of the minimum spanning tree. For each trajectory, coordinates were normalized by the maximal value for comparison with the other method.
Comparison of expression profiles among mouse lemur cell types
UMAP of cell types
To visualize similarities among the mouse lemur cell-type expression profiles, we embedded the high-dimensional 10x scRNA-seq expression data (around 30,000 genes) into a 2D UMAP. Cell types that were low quality (labelled with LQ in free annotation) or represented by fewer than 4 individual cells were excluded, which resulted in a comparison of 681 molecular cell types. Cell types were treated as pseudo bulk, with expression levels calculated for each gene by averaging the expression level of all cells within the same cell type and then taking the natural log transform (ln(avg count per 10K UMIs + 1)). Expression levels were further normalized by the maximal value of each gene across all cell types so that all ranged from 0 to 1. The cell-type gene expression matrix was then projected onto a 2D space with cosine distances between pairs of cell types used in the UMAP function (https://www.mathworks.com/matlabcentral/fileexchange/71902). Wilcoxon rank-sum tests were used to identify DEGs that distinguished related molecular cell types identified in the cell-type UMAP (for example, mature and progenitor sperm cells plus progenitor and proliferating immune cells versus proliferating cells of other compartments) as described in the Supplementary Methods.
Heatmap of cell-type pairwise correlation scores
To compare the overall gene expression profiles of cell types, Pearson’s correlation scores were calculated for every pair of cell types. Given data were obtained from different sequencing platforms (10x and SS2), we used the FIRM-integrated dataset as described above, which contains FIRM-generated PC coefficients for each cell. Cell types were treated as pseudo-bulk, and the cell-type average PC coefficients were calculated and used to determine the correlation coefficients. The cell-type pairwise correlation scores were plotted as a heatmap matrix (Extended Data Fig. 5d). Interactive forms of the heatmap matrix are available online at the Tabula Microcebus portal.
Evolutionary comparison of mouse lemur, human, macaque and mouse transcriptional profiles
Compiling comparable datasets
For the cross-species comparisons, we used published human, mouse and macaque scRNA-seq (and not single nucleus) datasets that were obtained using methods similar to those described above for the mouse lemur. These datasets included lung and muscle cells (from all compartments), epithelial cells of the liver, immune cells of the bone marrow and spleen, and germ cells of the testes (Supplementary Table 4). We manually re-annotated cells where necessary for consistency with the lemur annotations (see below). All lemur data were from the 10x data of this study with additional muscle data from L5 (ref. 22). Human data were from the 10x data of the Tabula Sapiens29, except for the lung, for which we used the 10x data from the Human Lung Cell Atlas14 and the testes, for which we used previously published drop-seq data35. Mouse data were from the 10x data of the Tabula Muris Senis13, except for the testes, for which we used previously published 10x data36. Given the limited data availability (either lack of the tissue or relevant cell types), we analysed only the lung and testes for the macaque data. Crab-eating macaque lung data were from previously published 10x data37 and rhesus macaque testis data were from previously published drop-seq data35. All datasets profiled adult animals; we excluded mouse postnatal developmental data from the analyses for consistency.
Orthology mapping across species
For orthology mapping, we merged the orthology databases from both NCBI and Ensembl (Supplementary Tables 5). We began by compiling all mouse lemur genes annotated in the NCBI (mouse lemur taxonomy ID: 30608), then merged the corresponding human and mouse orthologues from NCBI (gene_info.gz and gene_orthologs.gz from https://ftp.ncbi.nlm.nih.gov/gene/DATA/, February 2020). We next added Ensembl gene identifier (ID) numbers, gene names and human and mouse orthologue assignments from Ensembl Biomart (Ensembl Genes v.99, February 2020) using the Ensembl gene ID (variable ‘Gene_stable_ID’) for each NCBI gene ID (variable ‘NCBI_gene_ID’) in Ensembl Biomart. Mouse lemur genes that did not have an assigned human or mouse orthologue in either Ensembl or NCBI were removed, as were mouse lemur genes that had more than one human or mouse orthologue assigned, or that shared the same human or mouse orthologue with another mouse lemur gene. Note that unlike NCBI, Ensembl specifies the type of orthologue assignment (for example, ‘ortholog_one2one’ or ‘ortholog_one2many’); however, we did not use the Ensembl specification to filter one-to-one-to-one orthologues because, occasionally, a mouse lemur gene name was assigned by homology to multiple currently unnamed loci in Ensembl and because of this imperfect genome annotation, was labelled as sharing an ‘ortholog_one2many’ with human or mouse instead of ‘ortholog_one2one’. Finally, we appended the one-to-one orthologues between human and rhesus macaque and between human and crab-eating macaque, as assigned by Ensembl. A total of around 15,000 one-to-one gene orthologues were therefore uncovered across human, lemur and mouse genomes, around 14,000 across human, lemur, mouse and rhesus macaque genomes, and around 13,000 across human, lemur, mouse and crab-eating macaque genomes (Supplementary Table 5). Sequence identity was based on those reported in the Ensembl homology database.
Integrating cross-species datasets and unifying cell-type annotations
For the cross-species comparisons, we used the one-to-one gene orthologues that existed in all relevant datasets. Orthology mapping for the datasets was based on the NCBI or Ensembl gene ID if the original datasets provided the respective gene ID, and on the gene symbol if the gene ID was not provided. The choice of NCBI versus Ensembl depended on which version of the genome annotations the original dataset was aligned to. Some of the one-to-one orthologues were missing from one or more of the datasets; therefore, these were removed from the cross-species comparison. Together, we identified around 13,000 genes for the comparisons across human, lemur, and mouse genomes, and around 12,000 genes for the comparisons that also include either of the macaque species.
To unify cell-type annotations, human, mouse and macaque datasets were first re-annotated separately for each tissue and species using the same pipeline and marker genes as for the lemur data. For the male germ cells that formed a molecular gradient, we simplified the annotations into three discrete stages (spermatogonia, spermatocytes and spermatids) based on their original annotations and applied trajectory analysis (see below). Next, to ensure consistency of cell annotations across species, we applied Portal38 to integrate data from different species. Through adversarial learning of neural networks, Portal projects data into a space that minimizes species differences, from which an integrated UMAP is generated to visualize cell clustering from different species. Portal integration was performed separately for each tissue, except for bone marrow and spleen, which were jointly integrated. We manually inspected each integration UMAP and ensured that cells of the same designation showed reasonable cross-species co-clustering and separation from other cell types. We also made minor modifications to the cell annotations during this process to unify designations across species. For example, proliferating cells might co-cluster with the main non-proliferating population of the same cell type in the original dataset if the number of proliferating cells were too few (and they thus could not be distinguished by separate annotations), but they often formed a separate cluster with the proliferating populations of the other species in the integrated UMAP. In such a scenario, we re-annotated these cells as a proliferating subtype. We also merged cell types that had unclear cross-species correspondence and were almost indistinguishable in the species-integrated UMAP (for example, proliferating T, NK and NKT cells were grouped together and designated NK/T cells (PF)).
As additional validations of annotation consistency across species, we applied SAMap39,49, a self-assembling manifold algorithm and graph-based data integration method, to the lung and muscle datasets in order to identify orthologous (reciprocally connected) cell types on the basis of shared expression profiles across species. Cross-species cell-type similarity (visualized by the edge width in Extended Data Fig. 6d,e) is defined as the average number of cross-species neighbours of each cell relative to the maximum possible number of neighbours in the combined manifold. The default SAMap parameters were used in the analysis, and similarity scores less than 0.1 were removed.
Identifying species-unified trajectories
Trajectories were calculated for spermatogenesis across human, macaque, lemur and mouse datasets, as well as for three myeloid lineages (neutrophil, monocyte–macrophage and erythroid) of haematopoiesis across human, lemur and mouse datasets (macaque data not available). Trajectory detection and cell alignment was performed using the same custom in-house program as described above (in trajectory analysis), with the species-integrated UMAPs as input.
Calculating cross-species similarity scores for each cell type
Cell types with more than 15 cells in each of the species were used for the cross-species comparisons. This resulted in a total of 63 orthologous cell types for the comparisons across human, lemur and mouse data (63 × 3 = 189 total cell-type entries across all species), and 18 cell types for the comparisons of the lung cell types across human, macaque, lemur and mouse data (18 × 4 = 72 cell-type entries). Cell-type mean gene expression was calculated for each gene. Single-cell expression levels in the species-integrated dataset were normalized and log-transformed the same way as described above for the lemur-only dataset. That is, ln(UMIg/UMItotal × 1 × 104 + 1). Note, however, that because there were fewer genes in the cross-species dataset (only one-to-one orthologues), the absolute expression levels were higher than that in the lemur-only dataset.
We used correlation coefficients as a proxy for cross-species similarity. To score similarity for individual cell types, we calculated Spearman rank-based correlation coefficients of cell-type mean expression levels between human and lemur, rcHL and between human and mouse, rcHM. The cell-type mean expression levels were thresholded at ≥0.4 to mitigate the effect of background noise. Cross-species similarity was similarly calculated for cells at different stages and lineages of spermatogenesis and haematopoiesis by applying a moving window along the respective trajectories.
Calculating cross-species similarity scores for each gene
Cross-species similarity was calculated separately for individual genes using the tissue and three species (human, lemur and mouse) integrated dataset across the 63 orthologous cell types. We first quantified the mean expression (Emax) in the maximally expressed cell type in each species. Next, we filtered genes that were not expressed or expressed at low levels across the analysed cell types, requiring Emax > 0.5 in all three species, or Emax > 0.1 in all three species with Emax > 1.5 in at least one species. This resulted in a total of 7,787 genes for follow-up analysis. Mean cell-type expression levels across the 63 cell types were then normalized by Emax in each species and Pearson’s correlation coefficients between human and lemur (rgHL), human and mouse (rgHM) and lemur and mouse (rgLM) were calculated. We then calculated Δrg = rgHL – rgHM for each gene and tested the P value for Δrg being significantly higher (right-tailed) or lower (left-tailed) than 0 (see below). To identify genes with human–lemur-conserved but mouse-divergent expression patterns (that is, HL genes), we applied a threshold of Δrg > 0.4 and a right-tailed P value < 0.05. We also identified human–mouse-conserved and lemur-divergent (HM) genes and lemur–mouse-conserved and human-divergent (LM) genes using the same threshold levels. A similar analysis was also performed to detect genes that showed species-conserved or species-diverged expression patterns along the spermatogenesis trajectory and the neutrophil lineage of the haematopoiesis trajectory. We also identified genes that are highly conserved (rg > 0.8), lowly conserved (rg < 0.3) or moderately conserved (rg > 0.3 and rg < 0.8) in all three species. A full list of analysed genes and their statistics is provided in Supplementary Table 7. Expression patterns of example genes is visualized in Supplementary Fig. 4. Gene set enrichment analysis was performed using gProfiler50 for the highly conserved genes and HL genes, with all the analysed genes provided as a custom background gene set and otherwise default parameters.
To test whether one correlation coefficient (r) was significantly higher or lower than the other, we estimated the significance of their difference (Δr) being larger or smaller than 0 through Fisher’s Z-transformation. In essence, correlation coefficients, which were bounded and not normally distributed, were Fisher’s Z-transformed to the unbounded and approximately normally distributed space using the inverse hyperbolic tangent function, and their difference and respective P value were calculated using standard one-tailed t-tests in this transformed space. For display purposes, the mean and 95% confidence intervals were then inverse transformed and displayed in Fig. 5a and Extended Data Figs. 9a and 10a. Note that this inverse transformed Δr, which is bounded between −1 and 1, does not necessarily equal the initial Δr, which is between −2 and 2.
Identification of genes with primate-selective expression for each cell type
Using the cross-species dataset across human, lemur and mouse, we performed two separate Wilcoxon rank-sum tests for each gene and for each of the 63 orthologous cell types. The first was a two-tailed test comparing expression in lemur versus human, lemur versus mouse and human versus mouse. The second was a one-tailed test comparing expression in a cell type versus the rest of the cell types in the dataset (independently of the species). We calculated the fold change in mean expression for the above comparisons. Next, for each cell type, we searched for three categories of genes. First, genes with significantly primate-enriched expression, which requires that (1) cell-type mean expression of the gene is above 0.5 in both humans and lemurs and (2) 5-fold greater expression and P < 1 × 10–5 in both species compared with the orthologous mouse cell type. Second, genes with significantly primate-depleted expression, which requires that (1) cell-type mean expression of the gene is above 0.5 in mouse and (2) 5-fold lower expression and P < 1 × 10–5 in human and lemur compared with the orthologous mouse cell type. Third, genes that are significantly enriched in a cell type, which requires that (1) cell-type mean expression of the gene is above 0.5 in all three species and (2) 5-fold greater expression and P < 1 × 10–5 when comparing this cell type versus other cell types. The full list of the identified genes is provided in Supplementary Table 6.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
