key: cord-0264120-5l9psncw
authors: Liu, Shixuan; Ezran, Camille; Wang, Michael F. Z.; Li, Zhengda; Long, Jonathon Z.; De Vlaminck, Iwijn; Wang, Sheng; Kuo, Christin; Epelbaum, Jacques; Terrien, Jeremy; Krasnow, Mark A.; Ferrell, James E.
title: An organism-wide atlas of hormonal signaling based on the mouse lemur single-cell transcriptome
date: 2021-12-14
journal: bioRxiv
DOI: 10.1101/2021.12.13.472243
sha: 2c3475958d486a3a61f328034f6d27c4376ff430
doc_id: 264120
cord_uid: 5l9psncw

Hormones coordinate long-range cell-cell communication in multicellular organisms and play vital roles in normal physiology, metabolism, and health. Using the newly-completed organism-wide single cell transcriptional atlas of a non-human primate, the mouse lemur (Microcebus murinus), we have systematically identified hormone-producing and -target cells for 87 classes of hormones, and have created a browsable atlas for hormone signaling that reveals previously unreported sites of hormone regulation and species-specific rewiring. Hormone ligands and receptors exhibited cell-type-dependent, stereotypical expression patterns, and their transcriptional profiles faithfully classified the discrete cell types defined by the full transcriptome, despite their comprising less than 1% of the transcriptome. Although individual cell types generally exhibited the same characteristic patterns of hormonal gene expression, a number of examples of similar or seemingly-identical cell types (e.g., endothelial cells of the lung versus of other organs) displayed different hormonal gene expression patterns. By linking ligand-expressing cells to the cells expressing the corresponding receptor, we constructed an organism-wide map of the hormonal cell-cell communication network. The hormonal cell-cell network was remarkably densely and robustly connected, and included classical hierarchical circuits (e.g. pituitary → peripheral endocrine gland → diverse cell types) as well as examples of highly distributed control. The network also included both well-known examples of feedback loops and a long list of potential novel feedback circuits. This primate hormone atlas provides a powerful resource to facilitate discovery of regulation on an organism-wide scale and at single-cell resolution, complementing the single-site-focused strategy of classical endocrine studies. The network nature of hormone regulation and the principles discovered here further emphasize the importance of a systems approach to understanding hormone regulation.

Hormones are chemical messengers that circulate through the bloodstream and control and coordinate functions of both nearby and distant cells. Since the discovery of secretin over a hundred years ago (W. M. Bayliss, 1902) , about 100 mammalian hormones have been identified (Takei et al., 2015) . Over the decades, these hormones have been purified, sequenced, and synthesized, their receptors have been identified, and the intracellular signaling pathways they regulate have been explored (Tata, 2005) . Hormones regulate diverse biological processes ranging from growth and development, reproduction, metabolism, and immune responses, to mood and behavior. Many hormones are also therapeutically significant. For example, the discovery of insulin (Lewis and Brubaker, 2021) changed type I diabetes from an inevitably-fatal disease to a manageable one. Likewise, corticosteroids, a class of small molecule adrenal hormones, are used to treat diverse pathological conditions, including the inflammation attendant to COVID-19, where corticosteroids reduce the death rate in severe cases by ~30% (Horby et al., 2021; Sterne et al., 2020) .

The classical endocrine hormones are produced by secretory cells residing in one of the nine endocrine glands. However, it is now clear that other cell types, like fat cells and leukocytes, secrete global regulators into the circulatory system (Scheja and Heeren, 2019). In addition, some endocrine functions, like menstrual cycles in human females, depend upon the interplay and feedback between multiple hormones and multiple cell types (Jabbour et al., 2006) . For these reasons, a complete understanding of endocrine physiology requires comprehensive, global approaches. Variations in hormone regulation and function in different species further complicate the situation. Although comparative studies have shown that many hormones are highly conserved across vertebrates, the hormone-producing and -receiving cells are sometimes rewired for different functions (Tata, 2005) .

The grey mouse lemur (Microcebus murinus) is a small non-human primate with a lifespan of 5-10 years. Because of its small size, rapid reproduction, and close genetic proximity to humans, the mouse lemur is an appealing model organism for studies of primate biology, behavior, aging, and disease (Bons et al., 2006; Ezran et al., 2017; Kraska et al., 2011; Languille et al., 2012; Pifferi et al., 2019) . With the advances in single-cell RNA-sequencing (scRNAseq) techniques, we have recently created a molecular cell atlas of the mouse lemur by profiling transcripts from ~225,000 cells from 27 mouse lemur organs and tissues of 4 individuals (The Tabula Microcebus Consortium, 2021) (https://tabulamicrocebus.ds.czbiohub.org/). This organism-wide dataset provides an unprecedented opportunity to study global hormone regulation at single cell resolution in a non-human primate.

In this study, we present a comprehensive analysis of primate hormone signaling based on the information in the mouse lemur molecular cell atlas (The Tabula Microcebus Consortium, 2021) . We systematically mapped hormone-producing and target cells for 87 hormone classes, and created a browsable atlas that depicts gene expression of hormone ligands, modulators, and receptors. This systematic analysis detected both canonical and unreported sites of hormone secretion and action, and identified several examples of species-specific rewiring. We also analyzed the global expression patterns of hormone ligands and receptors, constructed and characterized an organism-wide cell communication network for hormone signaling, and systematically searched the network for feedback circuits. This analysis revealed a few remarkable principles of hormone signaling and a long list of potentially novel feedback mechanisms.

The hormone atlas provides a useful resource and a general approach to study hormone regulation on an organism-wide scale and at single-cell resolution, complementing the more reductionistic strategies typical in classical endocrine studies. The network nature of hormone regulation and principles Previous research has characterized tissue-specific functions of adiponectin in the liver, pancreas, skeletal muscle, endothelial cells, macrophages, and adipocytes (Ruan and Dong, 2016) . While we indeed observed expression of the adiponectin receptor ADIPOR1 in these tissues in the mouse lemur, we also found ADIPOR1 to be abundantly expressed in almost all lemur cell types (Fig. 1) . The ubiquitous expression of this receptor suggests that adiponectin exerts an organism-wide effect through ADIPOR1.

In comparison, ADIPOR2, the other adiponectin receptor, was expressed more selectively, most notably in the male germ cells and adipocytes. Interestingly, the mouse lemur expression pattern of ADIPOR2 differs from that of the mouse in which ADIPOR2 is reported to be most highly expressed in the liver and almost undetectable in the testis (Yamauchi et al., 2003) , suggesting differential ADIPOR2 effects in the two species.

Comparing expression patterns of the mouse lemur to that of the mouse and human, we identified interesting cases of hormone network rewiring. A telling example is motilin, canonically considered to be a gastrointestinal hormone secreted by enteroendocrine cells to regulate gastrointestinal smooth muscle movement (Kastin, 2013) . Mice lack both the motilin and the motilin receptor genes ( Table 2) . Indeed, genome comparative studies reported that motilin signaling is conserved in vertebrates but is lost in rodents (He et al., 2010) . The lack of motilin signaling in mice and rats has hampered progress on further characterizing motilin function as well as the development of motilin agonists for the treatment of gastroparesis (Sanger et al., 2013) . In the mouse lemur, we observed motilin expression in enteroendocrine cells in the small intestine but not colon (Fig. 1A) , which is consistent with the canonical view. However, there was also a broad range of other motilin-expressing cell types, including kidney nephron tubule cells, chemosensory tuft cells from both the airway (lung Prohormone processing enzymes. Prohormone processing enzymes are required for the maturation of many peptide hormones, including insulin, adrenocorticotropic hormone (ACTH), gastric inhibitory polypeptide (GIP), and somatostatin (Kastin, 2013) . The required processing enzymes have been identified for some but not all hormones (Seidah, 2011; Zhang et al., 2010) . In this study, we considered both the ligand and the known or suggested processing enzymes when determining if a cell type was a hormone source. Selective expression of the processing enzymes can help distinguish hormones that share the same precursor gene. For example, proglucagon (GCG) is the precursor for both glucagon and glucagon-like peptide-1 (GLP-1). GCG was expressed in both pancreatic α cells and intestinal enteroendocrine cells, but the two cell types expressed different processing enzymes: α cells expressed only PCSK2 which is required for glucagon synthesis, and enteroendocrine cells expressed only PCSK1 which is required for GLP-1 synthesis ( Fig. S1-glucagon) . Thus, as is the case in mice and humans (Sandoval and D'Alessio, 2015) , α cells are sources of glucagon and enteroendocrine cells are sources of GLP-1 in the mouse lemur.

Interestingly, we noticed that some cell types in the mouse lemur transcribed prohormone genes without expressing the respective processing enzymes. For example, proopiomelanocortin (POMC) is the precursor of multiple hormones including ACTH, melanocyte-stimulating hormone (MSH), βendorphin, and lipotropin, and the first step in POMC processing is carried out by PCSK1 (Wardlaw, 2011) . As expected, ACTH-secreting corticotrophs expressed both POMC and the necessary processing enzymes, but mouse lemur male germ cells expressed abundant POMC without PCSK1 ( Fig. S1 proopiomelanocortin-derived hormones). Similarly, male germ cells expressed GIP, which also requires PCSK1 processing. We also found renin (REN) expression in the intestinal and airway tuft cells without the necessary PCSK1, oxytocin (OXT) expression in kidney proximal tubule epithelium without the required carboxypeptidases E (CPE), and neuropeptide Y (NPY) expression in bladder urothelial cells without the necessary CPE (Fig. S1 ). These results suggest several possible hypotheses: different peptidases are involved, the prohormone transcripts are not translated, the prohormone proteins may function without proteolysis, or the processing enzymes may be expressed below detection levels in the animals that were sampled.

In summary, we have depicted the organism-wide, single-cell expression patterns for hormones, synthases, hormone modulators, and receptors from 87 hormone groups, and have systematically identified their sources and targets (Fig. S1) . This global analysis detected canonical cell types underlying hormone secretion and target functions while also revealing previously unreported sources and targets. This organism-wide single-cell hormone atlas should facilitate the systematic study of 8 hormone regulation in an unbiased, global manner, paving the way for the discovery of novel aspects of hormone regulation and function.

Next, we clustered mouse lemur cell types based on the transcriptional profiles of only the hormonerelated genes ( Fig. 2A-B) . We organized the cell types into 75 clusters and summarized the hormones and hormone receptors that were expressed in the majority of the cell types for each cluster in Table 3 . We noted that cells of the same type almost always clustered closely irrespective of their tissue of origin, and different cell types were well-separated in the cluster dendrogram ( Fig. 2A) . For example, all lymphatic endothelial cells clustered closely and showed characteristic IGF1 expression despite their coming from nine different tissues and four individuals (cluster 53, Fig. 3B, Fig. S1 -insulin-like growth factor, Table 3 ). Schwann cells from nine tissues clustered into myelinating and non-myelinating subtypes (cluster 40 and 37). Adipocytes sampled from four adipose depots and six other tissues formed one cluster (cluster 13, Table 3 ). Skeletal muscle stem cells from the limb muscle and tongue clustered together (cluster 17, Table 3 ). Mesothelial cells sampled from eight visceral organs formed a cluster distant from other stromal cells and showed distinctive expression of NPR1, CFB, PTGIS, ITGB3, and PNMT (cluster 44, Fig. 2A-B, Fig. S4A , Table 3 ). Although the analysis only examined the transcripts of 329 hormone-related genes, ~1% of all detected genes, these results suggest that hormone-related gene expression is indicative of the overall transcriptome among different cell types.

To further test this hypothesis, three comparative analyses were performed. First, we projected the highdimensional transcriptional profiles of the mouse lemur cell types into two-dimensional uniform manifold approximation and projection (UMAP) representations, which illustrate cell type similarity in the space of the selected genes. We compared the spatial relationship of cell types in two UMAPs: one constructed based on only the 329 hormone-related genes (Fig. 2B) , the other based on all genes excluding the hormone genes (>30,000 genes) (Fig. S3) . In each case, the data points represent cell types rather than individual cells (Fig. 2B, Fig. S3, Fig. S4 ). The two UMAPs showed remarkable similarity and cell types were similarly clustered with only minor differences (see Fig. S4C for comparison). One such difference was the chemosensory tuft cells that reside in the epithelium of both the airway (lung, trachea) and small intestine. In the space of hormone-related genes, tuft cells from all three tissues formed a cluster (cluster 21) with distinctive expression of renin (REN), cyclooxygenase 1 (PTGS1), arachidonate 5-lipoxygenase-activating protein (ALOX5AP, required for leukotriene production), and motilin (MLN) (Fig. 2B, Fig. S1 , Table 3 ). However, in the space of all genes, the lung, trachea, and intestinal tuft cells were separated from each other (labeled as 21*A, 21*B, and 21*C in Fig. S3 ). Another notable exception was skin melanocytes (cluster 38), which clustered closely with oligodendrocytes (cluster 39) when using only the hormone-related genes ( Fig. 2B ) but clustered with skin epithelial cell types (cluster 36) in the space of all genes (Fig. S3) . These differences notwithstanding, the hormone-gene and transcriptome-based cell type UMAPs were qualitatively similar and shared an overall clustering pattern.

We then compared the transcriptional distances between pairs of cell types for either just the hormonerelated genes (hormone-based distances) or the non-hormonal transcriptome (transcriptome-based distances) (Fig. 2C, Fig. S5A-B) . The two distances were positively correlated, with a Pearson correlation coefficient (r) of 0.70 (p < 1E-15). Similarly, we compared hormone-based transcriptional distances between pairs of cell types to a Cell Ontology-based distance. Cell Ontology distances represent a canonical view of the relatedness of cell types (Smith et al., 2007) and were calculated as the mean of the Cell Ontology embedding-based distances from the Cell Ontology graph and the text-based 9 distances from the text descriptions of the Cell Ontology terms . These two distances were also significantly correlated, with r=0.45 (p < 1E-15) ( Fig. 2D, Fig. S5A and C) . Taken together, these analyses support the notion that the small number of hormone-related genes faithfully classifies most of the cell types of the mouse lemur.

We then examined if the expression of hormone ligands or receptors alone could distinguish cell types (Fig. S4) . Interestingly, for most cell type clusters, hormone ligands or receptors alone classified cell types similarly. This suggests that most cell types have stereotypical expression of proteins involved in both incoming and outgoing hormonal signaling. A small fraction of clusters were segregated from other clusters mainly by only the ligands or receptors. For example, clusters of pituitary neuroendocrine cells (cluster 7) and tuft cells (cluster 21) are distinguished mainly by hormone ligands, whereas clusters of the skeletal muscle stem cells (cluster 17) and lung endothelial cells (cluster 51) are mostly based on differences in hormone receptor expression.

In summary, we found that cells generally show stereotypical and cell-type-specific expression of hormone ligands and receptors. Thus, hormone secretion and sensing essentially labels the identity of the cells.

Although cells of the same type tended to have similar hormonal expression patterns, we identified a number of interesting exceptions where closely related cell types differed in their hormone signaling. Below we highlight several such examples in which we identified changes in hormonal signaling among cells at different development stages (spermatogenesis), cells from different tissue origins (vascular endothelial cells), cell types that are not discretely segregated but form a continuous gradient (kidney nephron and vasa recta), cells of different molecular subtypes (hepatocytes), and normal vs. pathological (tumor) cells (Fig. 3) .

Spermatogenesis. Spermatogenesis is a tightly regulated process that happens asynchronously and continuously in the male testis (Griswold, 2016) . Therefore, sampling the testis at one fixed time point can capture snapshots of the entire process. The atlas includes seven sperm and sperm progenitor cell types that are at different spermatogenesis stages as indicated by the expression of known stage markers (The Tabula Microcebus Consortium, 2021) . Interestingly, these cells exhibited spermatogenesis-stagedependent hormone signaling (Fig. 3A) . The seven cell types fall into three major stages: the undifferentiated spermatogonia, the spermatocytes (which undergo meiosis), and the haploid spermatids that elongate and form mature sperms. Comparing hormone-related gene expression among these three major stages (cluster 1, 2, 8), over 50 hormone ligands, enzymes, and receptors were differentially expressed (Fig. 3A) , suggesting an extensive shift in hormone signaling. Previous studies on spermatogenesis have mainly focused on testosterone and on the supporting cells rather than sperms and sperm progenitors (Ruwanpura et al., 2010; Smith and Walker, 2014) . The present data suggest that sperms and sperm progenitor cells are not only differentially regulated by a variety of hormones but may also secrete a changing array of signaling molecules as they mature. Some of these molecules have been recently suggested to be involved in sperm functions. For example, progesterone was known to activate a sperm-specific calcium channel prior fertilization to induce Ca 2+ influx, but the underlying progesterone receptor(s) remained unidentified (Lishko et al., 2011) . Two non-classical progesterone membrane receptors (PAQR5/9) were highly expressed in the lemur spermatids and therefore may mediate the process (Fig. 3A) . Interestingly, two additional membrane progesterone receptors (PAQR8, PGRMC1) were specifically expressed in the undifferentiated spermatogonia (Fig. 3A) , suggesting that they may contribute to regulation of early spermatogenesis by progesterone (Lue et al., 2013) . Insulin-like factor 6 (INSL6), a member of the secretory relaxin family, was selectively expressed in the lemur spermatocytes (Fig. 3A) . Supporting a function in these cells, knockout of INSL6 arrests spermatogenesis at late meiotic prophase in mice (Burnicka-Turek et al., 2009) . This analysis suggests that sperms and sperm progenitor cells may be directly engaged in hormone signaling with surrounding and potentially distant cells to affect spermatogenesis and sperm functions. Vascular endothelium. Endothelial cells, including artery, vein, capillary, and lymphatic cells, line the surface of the circulatory system and are all derived from a developmentally uniform lineage (Dyer and Patterson, 2010) . The mouse lemur cell atlas includes a large number of endothelial cells captured from different organs, providing an opportunity for cross-organ comparison. Although vascular endothelial cells in general shared similar transcriptional profiles ( Fig. 2A-B, Fig. S3) , the lung and brain vascular endothelial cells exhibited unique hormone signaling signatures (Fig. 3B) .

Endothelial cells of artery, vein, and capillary in the mouse lemur lung showed higher expression of hormone receptors VIPR1, CALCRL, and RAMP1 ( Fig. 3B ) compared to endothelial cells from other organs. VIPR1 is a receptor for vasoactive intestinal peptide (VIP), a strong bronchodilator and one of the most abundant signaling peptides in the lung (Groneberg et al., 2001) . Similar expression of VIPR1 is also observed in human but not mouse lung endothelial cells (Tabula Muris Consortium, 2020; Travaglini et al., 2020) . RAMP1 forms heterodimeric receptors with CALCR or CALCRL that are regulated by several hormone ligands including CGRP, adrenomedullin, amylin, and calcitonin, all of which have been reported to have vasodilatory effects (Takei et al., 2015) . Interestingly, RAMP1 is expressed in the lung endothelium of the mouse lemur but not in mice or humans, whereas RAMP2 is expressed in all three species and RAMP3 is only expressed in humans ( Fig. S1 -adrenomedullin) (Tabula Muris Consortium, 2020; Travaglini et al., 2020) . These results suggest that lung vascular endothelium is a unique hormonal target in comparison to vasculature elsewhere, possibly because of the distinctive nature of the pulmonary circulation (e.g., high volume, low pressure) compared to the systemic circulation. The lung endothelium also shows species variation in the hormone signature, perhaps reflecting the relatively recent appearance of lungs in evolution compared to other visceral organs.

In contrast, brain vascular endothelium was found to express higher levels of a few hormone ligands/synthases. Endothelin-3 (ET-3), a member of the vasoconstricting endothelin family, was selectively expressed in the brain vascular endothelium (Fig. 3B) . ET-3 peptide has been detected at various brain locations in rodents, but it remained unsettled what cells secreted the peptide and what physiological role it played (Barnes and Turner, 1997) . Our results suggest that ET-3 is synthesized by the endothelial cells in the brain, which could target the nearby EDNRB-expressing SLC7A10+ astrocytes and the EDNRA-expressing vascular smooth muscle cells and pericytes of the brain (Fig. S1endothelin) . Additionally, the brain endothelium also expressed high levels of PTGES and PTGS2/COX-2 ( Fig. 3B) , which are the synthases for prostaglandin E2 (PGE2). PGE2 has been suggested to regulate inflammation-induced fever in mice (Engström et al., 2012; Wilhelms et al., 2014) . Interestingly, brain endothelial cells also selectively lacked expression of IGFBP4, which was abundantly expressed in the endothelium of all other tissues (Fig. 3B) .

The expression of the detected hormone signature genes (RAMP1, CALCRL, and VIPR1 in lung, and EDN3, PTGES, and PTGS2 in the brain) was higher in endothelial cells than other cell types from the same organ (Fig. S1) . Thus, it is unlikely that these organ-specific endothelial expression patterns resulted from cross-contamination (signal spreading). Altogether, lung and brain endothelial cells exhibited organ-specific hormone signaling which likely regulates organ-specific physiological functions.

Kidney nephron and vasa recta. Although discrete cell types were common in the mouse lemur cell atlas, there were clear examples of spatially continuous cell types that displayed a continuum of gene expression. These were the kidney nephron epithelium and the vasa recta endothelium (The Tabula Microcebus Consortium, 2021) , and these cells also displayed gradients of hormone signaling along their spatial axes ( Fig. 3C-D) .

Nephrons are single-cell layered, one-directional tubules that serve as the structural and functional units of the kidney. We identified a large number of hormone-related genes that are specific to different sections of the nephron epithelium (Fig. 3C) . Some of these hormones target specific nephron sections to regulate urine formation. For example, vasopressin is well known for its antidiuretic effect, promoting water reabsorption in the collecting duct principal cells (Takei et al., 2015) . Indeed, the vasopressin receptor AVPR2 was expressed exclusively in the collecting duct principal cells. Interestingly, the collecting duct intercalated cells expressed another vasopressin receptor, AVPR1A (Fig. 3C) , which has been recently reported to mediate the regulation of acid homeostasis by vasopressin (Giesecke et al., 2019) . Although the kidney epithelium is not commonly considered to be an endocrine tissue, many hormone ligands/synthases were also expressed in different nephron sections, including angiotensin (AGT), adrenomedullin 2 (ADM2), motilin (MLN), and neuropeptide Y (NPY). This suggests that nephron epithelium may play multifaceted roles in hormone production in addition to urine formation.

We also identified interesting differences in hormone signaling in the vascular endothelium of different regions of the mouse lemur kidney. In particular, we identified a gradual shift in hormone signaling along the arteriole-to-venule-axis of the vasa recta endothelium (Fig. 3D) . The vasa recta align parallel to the nephron tubules in the kidney medulla and function collaboratively with the nephrons in urine formation (Pannabecker and Dantzler, 2006) . Fig. 3D shows the hormone-related genes that are significantly differentially expressed along the vasa recta axis. Notably, several genes involved in IGF signaling were differentially expressed: IGF1, which was expressed in the lymphatic endothelial cells but almost absent from other vascular endothelial types (Fig. 3B, Fig. S1 -IGF), was highly expressed at the venule end of the vasa recta; IGF2 and IGF1R exhibited the opposite pattern with higher expression in the arteriole end; and four IGFBPs (IGFBP3-6), which modulate IGF signaling, exhibited axisdependent expression patterns. This spatial organization of IGF signaling may be evolutionarily conserved, as a recent study on mouse kidney vasculature development shows a similar pattern of IGF1 expression in the mouse vasa recta (Barry et al., 2019) .

Hepatocyte subtypes. The mouse lemur cell atlas identified multiple novel molecular subtypes of known cell types (The Tabula Microcebus Consortium, 2021) . For example, the liver hepatocytes of the mouse lemur clustered into two discrete subtypes, APOB+ hepatocytes and PHYH+ hepatocytes. These two hepatocyte subtypes are likely independent of the known zonal heterogeneity, and a similar subtype division was observed in mouse hepatocytes (The Tabula Microcebus Consortium, 2021). Interestingly, a number of hormone-related genes were among the most highly differentially expressed genes in the two hepatocyte subtypes in both the mouse lemur and mouse (Fig. 3E ) (The Tabula Microcebus Consortium, 2021) . Hormones including glucagon, androgens, and prolactin appear to selectively target the APOB+ hepatocytes as indicated by differential expression of the corresponding receptors (GCGR, AR, PRLR). In addition, angiotensin (AGT), a canonical hepatic hormone, was expressed at a much higher level in the APOB+ hepatocytes, whereas no subtype difference was detected in the expression of hepcidin (HAMP) or IGF1, two other hormones abundantly produced by the liver. Hepatocytes are known to show sex-dependent hormone regulation (Waxman and O'Connor, 2006) ; however, the subtype differences in hormone gene expression were not sex-dependent (Fig. S6) .

Liver is a central organ of metabolism and plays important roles in hormonal regulation. The detected heterogeneity in hormone signaling suggests potentially important functional divergence between the 12 hepatocyte subtypes in both the mouse lemur and mouse. A deeper understanding of the mechanisms and regulation of hepatocyte hormone heterogeneity would be particularly relevant to the study of mouse lemur-specific physiology, like hibernation/torpor and seasonal body weight changes (Terrien et al., 2018) .

The atlas also contains cells that are pathologic, including cells from the lung tissue of one female mouse lemur that were hypothesized by gene signature analysis and confirmed by histology to be metastatic tumor cells of endometrial adenocarcinoma origin (The Tabula Microcebus Consortium, 2021). The tumor cells clustered with the epithelial cells of the uterus from a second female individual (also, as it happened, with confirmed endometrial adenocarcinoma), particularly the MUC16+ non-ciliated uterus epithelial cells, and shared characteristic expression of oxytocin receptor (OXTR) and inhibin beta B subunit (INHBB) (Fig. 3F, Fig. 2A-B, Fig. S3, Table 3 ). This suggests that MUC16+ uterus epithelial cells may have been the origin of the metastatic tumor in the lung. However, certain aspects of hormone signaling were likely to have been rewired in the metastatic tumor (Fig. 3F) . Notably, the tumor cells showed a loss of estrogen receptors (ESR1), a change associated with high grade and advanced stage endometrial cancers in humans (Backes et al., 2016) and more often reported in high mortality type 2 human endometrial cancers (Shen et al., 2017) .

Interestingly, the metastatic tumor cells also acquired new abilities in hormone production. For example, elevated expression of prostaglandin E2 synthase (PTGES2) was present in tumor cells (Fig. 3F ). PGE2 has been suggested to be a tumorigenic factor in many cancers (Nakanishi and Rosenberg, 2013) . A recent study identified that the PGE2 synthase PTGES2 was elevated in human endometrial cancer tissues and found that PGE2 promotes cell proliferation and invasion in in vitro cell assays (Ke et al., 2016) . Resistin was also found to be highly expressed in the tumor cells (Fig. 3F) . Resistin has been reported to promote invasiveness of breast cancer cells in vitro (Lee et al., 2016) . Thus, endometrial tumors may utilize autocrine resistin signaling to promote metastasis. The shift in hormone signaling identified in the lemur endometrial tumor is similar overall to that reported in the human. This suggests that the rewiring of hormone signaling may be an important step in tumor development and metastasis. Mouse lemurs may serve as a useful animal model to study human cancers that are not well modeled in mice such as endometrial cancers (see more in (The Tabula Microcebus Consortium, 2021)).

In summary, we have highlighted five examples where closely related cells differ in their hormone signaling. These closely related cells usually share the same developmental origin but diverge during cell differentiation or development. Cells likely undergo changes in hormone signaling during the process to interact with other cells and accommodate their specific functional requirements.

Long-range hormonal signaling connects both nearby and distant cells through the bloodstream and enables cross-organ communication. By linking the source cells that produce a hormone ligand to the target cells that express the corresponding receptor, we constructed a global cell communication network for mouse lemur hormone signaling (Fig. 4A) . Nodes of the network are cell types and edges are directed from hormone source to target cells.

The hormone network was remarkably densely connected; 51% of all possible directed edges were found to be present. In fact, in the graphical depiction of the network (Fig. 4A) , individual edges are hardly distinguishable because the edges are so dense. The density of the hormone network is comparable to that of the cytokine network (Frankenstein et al., 2006) , but much higher than the 13 densities of other characterized biological networks. For example, the network densities of the yeast, fly, and human protein-protein interaction networks are all estimated to be ~0.1% (Giot et al., 2003; Rolland et al., 2014; Yook et al., 2004) , and the yeast genetic interaction network density is ~3% (Costanzo et al., 2010) . The dense connection between cell types by hormonal regulation is also highly robust: 67% of all existing edges are connected by more than one ligand-receptor pair, and the network maintains high density connection when randomly removing 10% of the ligand-receptor pairs (Fig. S7A) .

The density of the hormone network is largely due to the fact that multiple hormone ligands and receptors are expressed in most cell types (Fig. 4B) . On average, cell types expressed 17.0 ± 7.9 (mean ± S.D.) receptors or receptor complexes. In comparison, fewer hormone ligands were expressed (2.8 ± 2.5, mean ± S.D.). Only a few cell types expressed 10 or more hormone ligands (e.g., enteroendocrine cells, pancreatic endocrine cells, certain neurons, and certain fibroblasts), and 38.3% of cell types expressed a single ligand or no ligand. This indicates that cells in general respond to more hormones than they produce.

The five neuroendocrine cell types of the anterior pituitary-corticotrophs, thyrotrophs, gonadotrophs, lactotrophs, and somatotrophs-play a vital role in endocrinology. These cells are components of regulatory cascades, and classically they are regarded as receiving regulatory signals from one or a few hormones produced in the hypothalamus, and in turn producing one or a few hormones to regulate a peripheral endocrine gland. Yet even these cells were found to express numerous receptors, ranging from 17 (corticotrophs) to 29 (gonadotrophs and somatotrophs) ( Table 4 ). These cells appeared to be regulated by growth factors, neuropeptides, peripheral and classical hypothalamic hormones, as well as neurotransmitters (e.g., dopamine) that were not on our hormone list. Likewise, they were found to express genes for several hormones in addition to the classical hormones (Table 4) . For example, corticotrophs expressed galanin (GAL), gastrin (GAS), and cortistatin (CORT), as well as the ACTH precursor (POMC) and its processing enzymes PCSK1, PCSK2, and CPE.

Taken together, this analysis suggests that virtually all cell types across the body are extensively engaged in cell communication via hormonal signaling and that the connections are likely evolutionarily robust to random mutations that disturb hormone regulation. The dense connection of the network resulted in a large number of highly connected nodes and no obvious network cores (Fig. 4A, Fig. S8 ). This differs from the commonly studied scale-free networks, which are dominated by a small number of highly-connected hub nodes (Barabási, 2003) . Studies on hormone regulation usually focus on one specific site of secretion and/or target; our analysis emphasizes the importance of a systems view in understanding hormone regulation.

We next examined properties of the hormone network edges. In contrast to the protein-protein interaction network and the genetic interaction network, the edges of the hormone communication network are directional and can be distinguished into different types according to the specific ligandreceptor pairing. The hormone network contained 332 different edge types that are selected pairings of 128 ligands and 166 receptors. A ligand or receptor may contain multiple genes as required for ligand synthesis/maturation and in receptor complexes (see Table 1 ).

Some ligands and receptors were expressed in highly specific cell types while others were more broadly expressed. We therefore quantified for each ligand or receptor a generality score to evaluate the range of its expression across different cell types in the atlas ( Fig. 4D-E, Fig. S9 , Table 5 ). Strikingly, a small 14 number of hormone receptors were globally expressed. For example, the top ranking receptor, adiponectin receptor 1 (ADIPOR1), was expressed in 93.4% of the cell types. Other broadly expressed hormone receptors included the renin receptor ATP6AP2 (90.4%), the membrane-associated progesterone receptor PGRMC1 (84.1%), the cortisol receptor NR3C1 (73.7%), the resistin receptor CAP1 (72.9%), the IGF1 receptor IGF1R (69.5%), the opioid growth factor receptor OGFR (65.9%), and the thyroid hormone receptor THRA in complex with RXRB (54.3%). While it has been recognized that some of these hormones exert organism-wide effects, our results suggest that they target not only a broad range of tissues but likely regulate a large portion or perhaps majority of the cells of the body.

In comparison, most hormone ligands and receptors were expressed in specific cell types; ~85% of the hormone ligands and ~60% of the hormone receptors were detected in less than 5% of the cell types. Certain hormones like follicle stimulating hormone (FSHB and CGA) were only detected in a single cell type. Only a few hormone ligands were relatively broadly expressed, including adrenomedullin (ADM, 18.6%), acylation stimulating protein (C3, 18.1%), IGF2 (17.9%), and resistin (RETN, 15.4%). The generality score followed a power-law distribution toward the tail where a small number of globally expressed receptors dominated the distribution ( Fig. S9B-C) . This resulted in a bimodal distribution of the network node outdegree (abundance of outgoing edges/downstream regulation) ( Fig. 4A and C, Fig.  S8 ). Cells expressing hormone ligand(s) to the high-generality receptors (e.g., adipocytes which secrete adiponectin) sent globally connected outgoing edges, and other cells that only expressed ligand(s) to specific receptors sent a much smaller number of outgoing edges. In contrast, the nodes' indegree (abundance of incoming edges/upstream regulation) was more homogeneously distributed because of the higher specificity of hormone ligand expression and large number of hormone receptors expressed in the cells ( Fig. 4A-C, Fig. S8 ).

Many hormones have more than one receptor, which may mediate different functions ( Table 1) . We therefore examined if receptors for the same ligand tend to be co-expressed or if their expression is mutually exclusive across different cell types. To quantitatively compare the target cell type profiles, we calculated a Jaccard index (JI, intersection over union) for each qualifying receptor pair, namely the ratio of cell types that co-expressed both receptors to the cell types that expressed either of the receptors. This analysis revealed that the majority of the receptor pairs had extremely low JIs: 25% (41/161) of the receptor pairs had a JI equaling zero, meaning no overlap in their target cell types, and 57% (92/161) of the receptor pairs had a JI of less than 5% (Fig. 4F) . A telling case is the adrenoceptors; although epinephrine (adrenaline) targets many cell types, different adrenoceptor subtypes were selectively employed in different cell types (Fig. 4F, Fig. S1 -epinephrine/norepinephrine). For example, adrenoceptors ADRA1D and ADRB3 shared no target cell types; ADRA1D was expressed broadly in vascular smooth muscle cells and pericytes of multiple organs, certain glial cells, and leptomeningeal cells, whereas ADRB3 mainly targeted adipocytes. Vasopressin receptors AVPR1A and AVPR2 also had less than 3% overlap in their target cell types; AVPR2 was predominantly expressed in the kidney collecting duct principal cells, while AVPR1A was highly expressed in the kidney collecting duct intercalated cells and broadly expressed in many stromal cell types ( Fig. S1-vasopressin) .

A few receptor pairs showed a tendency toward co-expression, suggesting possible interactions in cell signaling. Supporting the notion, one top ranked pair was the IGF receptors IGF1R and IGF2R (JI=0.40), which are known to interact in IGF2 signaling. IGF2R internalizes after IGF2 binding and clears IGF2 from the cell surface, and therefore functions as a competitive inhibitor of IGF2-IGF1R signaling (Ghosh et al., 2003; Iams and Lovly, 2015) . Two other top ranked pairs were NPY1R/NPY5R (JI=0.49), receptors for neuropeptide Y, peptide YY and pancreatic polypeptide (Larhammar and Salaneck, 2004) , and ADIPOR1/ADIPOR2 (0.39), receptors for adiponectin. Coexpression of NPY1R and NPY5R has been identified in rodent brain, and recent studies indicate that the two GPCRs may form heterodimers, resulting in modulated agonist binding properties (Gehlert et al., 2007; Kilpatrick et al., 2015) . Similarly, biochemical studies on ADIPOR1 and ADIPOR2 in cell culture systems have shown that the two receptors can form heterodimers promoting membrane expression of ADIPOR2 (Almabouada et al., 2013; Keshvari et al., 2013) .

Taken together, this analysis showed that the expression of distinct receptors for the same hormone tends to be mutually exclusive, suggesting functional divergence and specification among hormone receptor subtypes. The rare cases of co-expressed receptors appear to involve functional interaction or cross-regulation between the receptor pairs.

Hormone regulation often involves feedback control. We therefore searched the mouse lemur hormone cell communication network to systematically identify potential feedback regulation. For simplicity, we focused on two-node feedback circuits, namely cell type pairs that are connected by edges in both directions (Fig. 5A ). This systematic search revealed a vast number of feedback circuits; 29% of all node pairs formed feedback circuits, showing that feedback control is widespread in the hormone network ( Table 6) . We then tested if feedback circuits are enriched in the hormone network, i.e. more prevalent than would be expected by chance. We applied degree-preserving random rewiring (Maslov and Sneppen, 2002) to randomly permute edge connections without changing node degree (i.e., network edge number and distribution). Fig. 5B shows the distribution of the feedback circuit counts from such randomly permuted networks (n = 1000). The number of feedback circuits in the permuted networks was consistently lower than in the original hormone network. However, the difference was small (on average, only ~0.2%). Additionally, we randomly removed 10% of ligand-receptor pairs to test how perturbations to ligand-receptor binding affect feedback distribution of the network. This analysis revealed that the network feedback circuits remained abundant despite edge removal as cell connections usually involved multiple ligand-receptor pairs (Fig. S7B, Fig. 4B ). Taken together, this analysis suggests that feedback circuits are abundantly and robustly distributed as a result of the dense connection of the hormone network. Evolution may have further enriched feedback circuits in hormone regulation as suggested by the random rewiring analysis.

We next reviewed selected feedback circuits between endocrine cell types ( Fig. 5C-L) to compare with literature reports on hormonal feedback regulation. The analysis successfully detected well-known feedback regulation. For example, Fig. 5C shows the feedback loop between gonadotrophs in the anterior pituitary and the estrogen-producing granulosa cells of the ovary (fortuitously obtained from female perigonadal fat), a classical example of hormonal feedback control that is evolutionarily conserved across vertebrates (Dwyer and Quinton, 2019) . Consistent with this scheme, mouse lemur gonadotrophs expressed two highly-conserved gonadotropins, follicle stimulating hormone (FSH, a dimer of FSHB and CGA) and luteinizing hormone (LH, LHB/CGA dimer), which target granulosa cells through the corresponding receptors FSHR and LHCGR. Granulosa cells, in turn, expressed key synthases required for the production of sex hormones, which feed back to gonadotrophs through androgen and estrogen receptors (AR, ESR1). Additional feedback regulation from granulosa cells to gonadotrophs was also mediated through anti-Müllerian hormone (AMH) and inhibin/activin signaling. We also detected the well-established feedback loop between somatotrophs in the anterior pituitary and hepatocytes in the liver, which was mediated through growth hormone (GH), IGF1 (Martín-Estal et al., 16 2016; Takei et al., 2015) , and potentially hepcidin and FGF21 (Inagaki et al., 2008) signaling (Fig. 5D) . Additionally, in the pancreas, we detected the known feedback circuits between pancreatic α and β cells through glucagon and insulin signaling and between α and δ cells through glucagon and somatostatin signaling ( Fig. 5E-F) (Hartig and Cox, 2020) . The detection of these various known feedback loops supports the validity of the comprehensive scRNAseq approach.

The analysis also revealed a long list of potential endocrine feedback circuits that have not been reported, or have been reported recently but are not yet widely appreciated. For example, we detected a feedback loop between pancreatic α cells and hepatocytes with one leg from α cells to hepatocytes mediated by glucagon and the other leg from hepatocytes to α cells through IGF1 and FGF21 (Fig. 5G) . Regulation of hepatocyte glycogenolysis and gluconeogenesis by glucagon is well-established (Habegger et al., 2010) . Previous work on the fasting-induced hormone FGF21 has shown that injection of glucagon increases plasma FGF21 in humans (Habegger et al., 2013) and injection of FGF21 reduces glucagon in mice (Mu et al., 2012) . A recent study also observed that plasma IGF1 and glucagon levels are negatively correlated in non-diabetic humans and found that IGF1 negatively modulates glucagon secretion in an α cell culture model (Mancuso et al., 2017) . Considering that the liver is a major source of circulating IGF1 and FGF21, our results and the literature reports suggest that α cells likely form a negative feedback loop with hepatocytes in maintaining metabolic homeostasis.

Additionally, adipocytes, through adiponectin, formed feedback loops with pituitary neuroendocrine cells, including the growth hormone secreting somatotrophs, the prolactin secreting lactotrophs, and the ACTH secreting corticotrophs (Fig. 5H-J) . It has been established that growth hormone, prolactin, and ACTH reduce adiponectin secretion in adipocytes (Carré and Binart, 2014; Iwen et al., 2008; Kopchick et al., 2020) , but much less is known about the effect of adiponectin on pituitary neuroendocrine cells. Interestingly, a recent study using cultured primary pituitary cells from two non-human primate species found that adiponectin decreases growth hormone and ACTH release and increases prolactin release (Sarmento-Cabral et al., 2017) . This suggests that adipocytes form a negative feedback loop with lactotrophs which possibly stabilizes prolactin and adiponectin levels, and form double-negative feedback loops with somatotrophs and corticotrophs, possibly functioning as bistable switches (Ferrell, 2002) . It remains to be characterized how adiponectin affects pituitary hormone secretion, what are the physiological functions of the pituitary-adipocyte feedback control, and whether this regulation is specific to primates.

We also detected feedback circuits of intestinal enteroendocrine cells with pancreatic β cells and adipocytes ( Fig. 5K-L) . It is known that enteroendocrine cells secrete GIP, which induces glucosedependent insulin secretion in β cells after eating and exposure of the intestinal lumen to nutrients (Takei et al., 2015) . Despite the important physiological role of GIP in insulin secretion and glucose homeostasis, there is limited research on the potential feedback regulation of GIP secretion by insulin. García-Martínez et al. reported that insulin induces prolonged GIP expression in an enteroendocrine cell line (García-Martínez et al., 2014) , suggesting a positive feedback regulation and synergistic effects between the two hormones in maintaining glucose homeostasis under a surge of blood glucose after eating.

Taken together, this analysis identified an extensive list of potential feedback circuits in the global hormonal cell communication network (Table 6) , providing a rich resource for generating hypotheses on endocrine regulations and to compare hormonal regulatory mechanisms in different species.

Lastly, the mouse lemur hormone atlas may bring new insights to the study of mouse lemur physiology, such as circannual rhythms (Perret and Aujard, 2001a; Petter-Rousseaux, 1979) , social interactions (Perret, 1992) , and aging (Epelbaum and Terrien, 2020; Kraska et al., 2011; Languille et al., 2012; Pifferi et al., 2019) . As an adaptation to seasonal environmental changes, mouse lemurs have evolved striking seasonal rhythms in metabolism, behavior, and reproduction throughout their 5-10-year lifespan. Mouse lemurs build fat reserves in preparation for the winter and enter daily torpor, a hypometabolic state, to preserve energy, while in the summer the animals maintain a lean body and regrow their reproductive organs for mating. The global nature of these phenotypic changes suggests that they may be hormonally mediated. To date, only 12 circulating hormones have been measured in mouse lemurs, but all showed dramatic seasonal changes (Fig. 6A, Table 7 ). This suggests that hormones are likely the key messengers coordinating and driving various seasonal changes among different organs. The hormone atlas suggests that these seasonally changing hormones target a broad array of cells and tissues; 95% of the cell types possessed receptors for at least one of these 12 hormones. This finding suggests that seasonal rhythms engage not only tissues that demonstrate obvious morphological changes (e.g., adipose tissues and gonads), but involve global shifts in cellular metabolism and physiology (Fig. 6B ).

Hormonal signaling is vital to multicellular life, allowing cells to communicate with both surrounding cells and remote cells in other tissues. In this study, we systematically mapped hormone-producing andtarget cells across the 739 cell types from 27 different tissues from the mouse lemur molecular cell atlas. We generated a browsable atlas that depicts gene expression of ligands, modulators, and receptors for 87 hormone classes (Fig. S1, Fig. 1 ). This dataset complements the scale and resolution of classical endocrine studies and brings a cellular and molecular foundation for studying novel aspects of hormone regulation. The mouse lemur hormone atlas represents a primate model of hormone signaling and should be useful in future organism-wide cross-species comparisons.

The hormone atlas also provides an opportunity to extract global principles of hormone signaling. The hormone-related genes exhibited cell-type-dependent, stereotypical expression patterns, and their transcriptional profiles faithfully classified cell types, despite their constituting less than 1% of the transcriptome (Fig. 2, Fig. S3, Fig. S5 ). Ligands alone and receptors alone also clustered cell types well, with certain clusters mostly defined by their ligands and others by their receptors (Fig. S4) . It is remarkable, perhaps, how well a cell's intercellular communication repertoire encapsulates the cell's identity.

There were, however, some exceptions to this general rule. For example, sperms and sperm progenitor cells exhibited dramatic changes in hormone gene expression during spermatogenesis despite their close developmental connection (Fig. 3A) . In addition, vascular endothelial cells from lung and brain exhibited hormone signatures that were distinct from those of other vascular endothelial cells (Fig. 3B) . Hormone signaling varied in an orderly and graded fashion along the epithelium and vasa recta endothelium of the nephron ( Fig. 3C and D) . Hepatocytes were found to fall into two molecular subtypes with differential hormonal gene expression (Fig. 3E, Fig. S6 ), suggesting different physiological roles. And finally, metastatic endometrial tumor cells exhibited different hormonal signaling from that seen in their likely cells of origin in the uterus (Fig. 3F) , suggesting possible rewiring of hormone signaling during tumor development and metastasis.

The cell-type specific gene expression data allowed us to construct an organism-wide cell communication network for mouse lemur hormone signaling, by connecting the hormone-producing cells to the corresponding target cells (Fig. 4) . This hormone network was exceptionally dense and robust to deletion of individual edges (Fig. 4A, Fig. S7A ), which mainly resulted from two network features. First, most cells in the network expressed multiple hormone receptors and/or ligands, though in general cells were regulated by more hormones than they produced (Fig. 4B) . Second, although most hormone ligands and receptors were expressed in specific cell types, a small fraction of hormone receptors were nearly ubiquitously expressed (e.g., adiponectin receptors, a membrane progesterone receptor, the IGF1 receptor, the thyroid hormone receptor, and the cortisol receptor), resulting in almost global connection of the relevant hormone-producing cells to the other cells of the body (Fig. 4C-E, Fig.  S8, Fig. S9) . Interestingly, despite the density of the hormone signaling network, subtypes of hormone receptors that bind to the same ligand were generally expressed in a mutually exclusive fashion, supporting the notion that different hormone receptor subtypes usually mediate different physiological functions through separate target cells (Fig. 4F) . The dense connection of the network indicates extensive hormonal communication among the cells and contributes to a large number of highly connected nodes and no obvious network cores. Whereas regulatory circuits of individual hormones may be hierarchical, as classically viewed, the overall decentralized network architecture of the hormone interactome highlights the importance of taking a comprehensive approach to the understanding of global hormonal function.

Feedback regulation was widespread in hormonal signaling, and it was robust to network perturbations (Fig. 5B, Fig. S7B ). The analysis identified well-known feedback circuits among several endocrine cell pairs and detected additional, previously unreported, hormones involved in feedback control (Fig. 5C-F) . We also identified many novel feedback circuits (Fig. 5G-L) which may be involved in adaptation or homeostasis (negative feedback circuits) or in generating switch-like responses (positive and doublenegative feedback circuits). This hypothesis-generating analysis provided an extensive list of potential feedback circuits ( Table 6 ) to guide future functional validation and characterization studies. Our findings may also inform studies of the systems biology and network topology of hormone signaling (see for example (Karin et al., 2020) ).

Mouse lemurs display dramatic seasonal changes in body weight, hibernation propensity, and reproductive activities (Perret and Aujard, 2001a, 2001b; Petter-Rousseaux, 1979) . Such seasonal rhythms are likely regulated, at least in part, by seasonally-varying hormones that coordinate the changes across different organs of the body (Fig. 6A) . The hormone atlas shows that almost all (95%) of the mouse lemur cell types are regulated by seasonally varying hormones (Fig. 6B) . This suggests that most, if not all, mouse lemur cells across all organs of the body experience seasonal changes in metabolism, growth, and function. By systematic mapping of hormone secreting and target cells, the atlas also points to candidate cell types and signaling to explore the upstream and downstream regulation for each of the seasonally varying hormones. Circannual rhythms are widespread in nature (Gwinner, 2011; Lincoln, 2019) and humans display seasonally varying hormone levels and cellular transcriptomes (Dopico et al., 2015; Sailani et al., 2020; Tendler et al., 2021) . The findings of the present study may apply to species beyond the mouse lemur, and may bring new insights for season-associated human diseases (Kurlansik and Ibay, 2012; Marti-Soler et al., 2014; Moriyama et al., 2020) because of the close genetic proximity of mouse lemurs and humans.

Single-cell RNA sequencing is a powerful tool for assaying cellular function and activity, but there are limitations. Although the technique has improved rapidly in sensitivity, accuracy, and efficiency, scRNAseq still faces low capture efficiency and high transcript dropout . With the large number of cells assayed in the mouse lemur molecular cell atlas, we averaged gene expression across cells of the same cell type to smooth the noise in individual cells from dropout. In the analysis of hormone regulation, genes involved in hormone production, such as ligands and modulators, tend to be abundantly expressed; but certain hormone receptor proteins may be low in copy number and slow in turnover, and therefore may result in false-negative errors. In addition, there can be false-positive errors caused by contamination from other cells in the same sample (i.e., signal spreading). This is most notable for abundantly expressed ligand genes and tissues affected by endogenous digestive enzymes (e.g., pancreas). In this study, we scored positive expression of a hormone ligand only when both the ligand gene and the necessary modulators (such as prohormone processing enzymes) were present. This approach mitigated false-positive detection of highly abundant ligands. Lastly, scRNAseq has low capture efficiency for certain cells such as neurons and large sized adipocytes. Nevertheless, the mouse lemur cell atlas included several neuronal subtypes from the brain cortex, brainstem, and retina, as well as both UCP1-high and UCP1-low adipocytes from multiple tissues (The Tabula Microcebus Consortium, 2021).

Although not yet a complete sampling, the mouse lemur cell atlas is one of the most comprehensive and accurately-annotated organism-wide scRNAseq datasets reported to date. Other organism-wide scRNAseq atlases have been created or are being created for other species, including mouse (Han et al., 2018; Tabula Muris Consortium, 2020 ; The Tabula Muris Consortium, 2018), ascidian (Cao et al., 2019) , cynomolgus monkey (Han et al., 2020) , and human (The Tabula Sapiens Consortium, 2021), bringing exciting resources to the study of the evolution of hormone signaling, as well as the cellular and molecular mechanisms and interactions underlying hormone-regulated biology.

Traditionally, biologists focus on one particular cell type and/or organ system. However, hormone signaling by its very nature connects different cell types and organs. The emerging organism-wide atlases open a new door in understanding the cross-organ communication systems that allow organisms to function as a reliable robust whole.

Mapping and counting transcript reads: The single-cell gene expression matrix of the mouse lemur cell atlas was obtained from (The Tabula Microcebus Consortium, 2021) where methods and parameters to extract the transcript counts can be found. In this study, we used only the 10x droplet-based sequencing derived data for all analysis. In brief, the Microcebus murinus genome assembly (Mmur 3.0, NCBI) was used to create a genome FASTA file. Sequencing reads were processed with Cell Ranger (version 2.2, 10x Genomics) to generate cell barcodes and the gene expression matrix. We additionally counted reads for three genes (MC1R, SCT, LHB) that are unannotated in Mmur 3.0, NCBI. For MC1R and SCT, we used their annotation in the Ensembl Microcebus murinus genome. For LHB, we used the uTAR-scRNA-seq pipeline described by (Wang et al., 2021) to first predict its chromosome location in the mouse lemur genome. In brief, the pipeline employs the approach of (Chae et al., 2015) to detect transcriptionally active regions (TARs) from aligned sequencing data and annotates these TARs as unannotated (uTARs) or annotated (aTARs) based on their overlap with the existing annotation. For the LHB_uTAR counts, we uncovered uTARs from the mouse lemur pituitary and hypothalamus dataset and used BLASTn (Altschul et al., 1990; Johnson et al., 2008) to find nucleotide sequence homology of all uTARs against human LHB. We found one uTAR at position NC_033681.1:19143249-19145949 with high sequence homology to human and macaque LHB. The LHB_uTAR is highly differentially expressed in the expected cell type (gonadotroph). The LHB_uTAR is located next to gene RUVBL2, in conservation with the co-location pattern of the two genes in humans and mice (Parfait et al., 2000) . We therefore conclude that the identified LHB_uTAR is the likely coding region of the mouse lemur LHB 20 gene. We then used Drop-seq tools (Macosko et al., 2015) to extract transcript counts for the three genes (MC1R, SCT, LHB_uTAR).

We compared the sequence similarity between orthologs of human and mouse or lemur. Protein sequences were retrieved from NCBI according to protein RefSeq using the Matlab built-in function getgenpept. Conversion between gene ID and protein Refseq identifier was performed using bioDBnet (https://biodbnet.abcc.ncifcrf.gov/db/db2db.php). We aligned ortholog protein sequences and used the maximal scores in the cases of one-to-multi or multi-to-multi ortholog mapping. Global alignment was carried out using the Needleman-Wunsch algorithm and default parameters of the function nwalign.

The functional concentration is different for different proteins, so we applied adaptive thresholding to determine if a gene is positively expressed in a cell type, where cell type is defined as a unique combination of the cells' annotation and tissue of origin, as described in (The Tabula Microcebus Consortium, 2021) . For example, vein cells of the lung and brain were considered two different cell types. Cell types that were of low quality (low gene count or transcript count) or annotated as technical doublets were removed from the analysis. For each gene, two cell type parameters were used, average expression level and percentage of nonzero cells. Single-cell expression levels were first normalized by total transcript counts, multiplied by 10,000, averaged across all cells in the cell type, then natural log transformed (i.e., log(1+AverageCountPer10KTranscripts); note that 1 is added to the average so that genes with zero transcripts will yield zero rather than negative infinity). To determine positive expression at the cell type level, we first applied a low-level absolute threshold on average expression level (0.037) and percent of nonzero cells (5%) below either of which we considered the expression to be negative. Second, we calculated a "maximal" expression level of the gene by searching the high-expression levels across the atlas. Specifically, we considered the 99th percentile of the cells that passed the first thresholding to be the maximum to increase the robustness of the analysis. Lastly, we applied relative thresholding of the "maximal" level (for ligand genes at 20% and other gene types at 5%). Cell types that passed all thresholding criteria were considered to have functional expression of the gene. We used a higher relative threshold for the ligand genes because 1) ligands are secreted into the circulation so a highly expressing source outshadows a low expressing source, whereas non-ligand products of other genes (e.g., receptors, synthases, processing enzymes) function inside individual cells; 2) we noted that ligand genes were usually expressed at higher levels and were more easily affected by cross-contamination (signal spreading) inside a tissue. For hormone ligands and receptors that included multiple genes, we applied AND or OR logic gates as indicated in Table 1 to determine if a cell type expressed the ligand or receptor. For an AND logic gate, a functional cell type should express all its components above the functional level; while for an OR logic gate, expression of any of its components would suffice.

We performed hierarchical clustering of the cell types included in the mouse lemur atlas based on the transcriptional profiles of hormone ligand, synthase, modulator, and receptor genes. Analysis was carried out at the cell type level and relative average gene expression levels were used to calculate cell-type pairwise distances. To test consistency of hormone-related gene expression among mouse lemur individuals, we separated cells from different lemurs. Because circulating immune cells sampled from different tissues were already kept as different cell types, immune cell types were not separated for different lemurs to avoid overcrowding of the data. 21 We found that this simplification only slightly influenced the clustering results of individual cell types and did not affect the overall clustering pattern. Next, cell types were filtered to remove potentially low quality cell types that were: 1) low in cell number, 2) technical doublets, or 3) low in transcript and/or gene count. We also manually corrected for notable cross-contamination in certain cell types. For example, almost all non-endocrine cell types of the pancreas showed high expression of insulin (INS) and proglucagon (GCG), probably caused by cross-contamination due to autolysis in the pancreas during tissue processing. We therefore removed INS and GCG expression from all non-endocrine pancreatic cell types. Similarly, we corrected contamination of GUCA2A, GUCA2B, GIP, NTS, and SST from stromal, endothelial, and immune cells of the intestine (small intestine and colon). Matlab built-in functions were used in the hierarchical clustering analysis. We used cosine distances to calculate cell type pairwise distances and then applied average-linkage to construct the dendrogram (cophenetic correlation coefficient = 0.80). Cell type clusters were then divided at a threshold of 0.49 with minor manual inputs to further divide certain large clusters or merge an apparently close neighboring leaf. We obtained 75 total clusters and numbered them according to their sequence in the dendrogram shown in Fig. 2A . We also named the clusters according to the shared characteristics of its cell types. Cell types included in each cluster are listed in Table 3 with the list of hormones, hormone receptors, and hormone-related genes that were expressed in the majority (≥ 50%) of the cell types for each cluster.

In addition to visualizing the cell type clustering via dendrogram, we also applied uniform manifold approximation and projection (UMAP) (McInnes et al., 2018) to embed the high dimensional expression data to 2D using a package developed by (Meehan et al., 2021) . Cell type pairwise cosine distances were used as input when generating the UMAPs.

For comparative analysis, we additionally calculated distances between cell types by either the nonhormonal transcriptome or by Cell Ontology (Smith et al., 2007) . Transcriptome-based distances were calculated similarly as the hormone-based distances were, with relative expression levels of nonhormonal genes used to calculate pairwise cosine distances. For Cell Ontology distances, each cell type was assigned a Cell Ontology term during cell type annotation as described in (The Tabula Microcebus Consortium, 2021). Cell Ontology distances were calculated as the mean of Cell Ontology embeddingbased distances and text-based distances . Cell Ontology embedding-based distances represent distances between cell types in the Cell Ontology graph. Text-based distances were calculated using an adapted version of Sentence-BERT (Reimers and Gurevych, 2019; Wang et al., 2020) .

To estimate the quality of the clustering results, we calculated silhouette values and compared cluster average silhouette values based on only the hormone ligand genes, hormone receptor genes, or the non- cell types belonging to a cluster with more than one member. a(i) is the mean in-cluster distance, namely the mean distance between cell type i and all other cell types in the same cluster. b(i) is the minimal outcluster distance, namely the minimal mean distance between cell type i and members of other clusters.

Vasa recta trajectory analysis: Analysis was performed independently for kidney samples from three mouse lemur individuals using an in-house program written in Matlab. To detect the arteriole-to-venule trajectory, data were first pre-processed following the standard procedure as in Seurat (Butler et al., 2018) . In brief, single cell transcript counts for all sequenced genes were normalized to total transcript counts, log transformed, and linearly scaled so that each gene has a mean expression of 0 and standard deviation of 1. Then genes with high dispersion (variance-to-mean ratio) compared to other genes with the same mean expression were selected as highly variable genes for principal component analysis (PCA). Following PCA, the top 20 principal components that were not driven by extreme outlier data 22 points or immediate early genes were used to generate a 2D UMAP using cell-cell Euclidean distances as input (McInnes et al., 2018) . To quantify the trajectory, we detected the probability density ridge of the data points in the UMAP using automated image processing. The direction of the trajectory was manually assigned according to marker gene expression. Individual cells were aligned to the trajectory by finding the shortest connecting point to the trajectory. Data points that were too distant from the trajectory were deemed outliers and removed from the following analysis. To detect the hormone-related genes that follow a trajectory dependent expression pattern, we calculated the Spearman correlation coefficient between the gene expression level and 20 unimodal patterns. These preassigned patterns have their data smoothly changing along the trajectory and have a single peak. We let the 20 patterns have their peaks uniformly distributed from 0 to 1. This analysis tested, for each gene, whether its distribution along the trajectory follows one or more of the 20 unimodal patterns. All p values were then multitesting corrected to detect significant genes and their associated patterns. Lastly, we compared results from the three kidney samples and kept the genes that showed consistent patterns in all three samples.

We used the Wilcoxon rank-sum test to detect differentially expressed hormone genes among cell populations. Genes that were expressed above functional level in any of the analyzed cell populations were included in the analysis. For hormonerelated genes differentially expressed during spermatogenesis, we tested whether each gene is differentially expressed in one of the three sperm cell clusters compared to the other two clusters. For the analysis of endothelial cells, we tested whether each gene is differentially expressed in the lung or brain endothelial cells compared to all other endothelial cells. For the kidney nephron analysis, we tested whether each gene is differentially expressed in one of the nephron clusters compared to the rest of the nephron epithelial cells. For the hepatocyte analysis, we tested whether each gene is differentially expressed in one subtype compared to the other in both the mouse lemur and the mouse. The mouse hepatocyte data was obtained from (Tabula Muris Consortium, 2020) smartseq2 plate-based sequencing of mouse hepatocytes. For the tumor cell analysis, we tested whether each gene is differentially expressed in the metastatic cells versus the MUC16+ non-ciliated uterus epithelial cells. All p values were multi-testing corrected, and then all genes with corrected p values less than 0.05 were examined to confirm consistent expression of the gene in all relevant cell types and mouse lemur individuals.

Cells with the same annotation name and tissue of origin were grouped as a cell type/node. Data from different lemur individuals were merged in the analysis. Cell types that were likely low quality were removed. Edge connections were determined independently for all possible cell type pairs. For example, if cell type A expressed a hormone ligand and cell type B expressed the corresponding hormone receptor(s), a directed edge from A to B was drawn. Note that this allows self-loops on any cell type that expresses both the ligand and corresponding receptor. When multiple directed edges connected two cell types, they were merged for the purpose of edge counting.

Network density is defined as the number of edges divided by the number of possible edges in the network. Because the hormone network is directed, its network density is calculated as 2 . For comparison, the density of an undirected network (e.g., protein interaction network) is calculated as

if it allows self-loops, or 2 * ⋅( −1) if it does not allow self-loops. 23 Outdegree and indegree are node features and defined as the relative number of outgoing edges (node expressing a ligand) and incoming edges (node expressing a receptor). Edge counts were normalized by cell type abundance so that node degrees range from 0-100% and represent the fraction of the cell types that the node connects with by outgoing or incoming edges. To normalize the node degrees, we used a weighted cell type abundance determined by the clustering so that each cluster contributes the same weight. This is because certain cell types like circulating immune cells are likely over-represented in the network as they appear in multiple tissues and would be categorized as different "cell types" by our definition. We did not merge circulating cell types from different tissues because it is challenging for some cell types like T cells to determine whether they are circulating or potentially tissue specific or tissue-resident. We also tried simply normalizing the node degree by the total number of nodes; similar conclusions were reached using either method.

To calculate generality scores for the hormone ligands and receptors, we counted the cell types that express the ligand or receptor. Cell types that were likely low quality were removed from the analysis. Cell type counts were normalized by cell type abundances similarly as was the node degree, so that the generality score was not dominated by over-represented cell types. In comparison, we also counted the cell type clusters that express the hormone ligands and receptors. Here we defined that a gene was expressed in a cell type cluster if more than 50% of the cell types in the cluster positively expressed the gene. For hormone ligands and receptor complexes that contained multiple genes, we applied the AND or OR logic gates indicated in Table 1 to determine positive expression in cell type clusters.

To analyze the distribution of generality scores, we tested if the distribution followed normal, exponential, lognormal, or gamma distributions. We first performed parametric distribution fitting and then performed Kolmogorov-Smirnov tests comparing the fitted versus data cumulative distribution functions. This analysis used Matlab built-in functions fitdist and kstest. Additionally, we tested if the data followed a power law distribution using plfit and plpva functions developed by (Clauset et al., 2009 ).

We calculated Jaccard indices for pairs of receptors that 1) can bind to the same ligand as described in Table  1 , and 2) were both expressed in at least one atlas cell type. Jaccard index here is calculated as the number of cell types that express both receptors divided by the number of cell types that express either of the receptors. Cell types that were likely low quality were removed from the analysis. Data from different mouse lemur individuals were merged in the analysis.

We performed two types of network perturbations: degree-preserving random rewiring and random edge removal. For degree-preserving random rewiring, we applied the algorithm developed by (Maslov and Sneppen, 2002) . The algorithm starts with the original network and performs repeated rewiring steps. For each rewiring, it randomly identifies two edges (e.g., A→B and C→D) and swaps the connections (i.e., A→D and C→B). This allows constructing a random network while preserving the node indegree and outdegree of the original network. For the perturbation by edge removal, we randomly selected 10% of the edge types (ligand-receptor pairs) and removed these edges from the original network. Note for this analysis, multiple directed edges connecting the same cell type 24 pairs were not merged. Therefore, the connection between two cell types would be removed only if all of its redundant edges (of different ligand-receptor pairs) were selected to be removed. Table 3 . The lower part of the panel shows the heatmap of the relative gene expression levels for the hormone-related genes. Cluster-specific hormone genes are ordered based on the cluster with the highest expression level, with the genes preferentially expressed in the left-hand cell type clusters on the top. Broadly expressed genes (expressed in more than 35% of cell types) are placed in the bottom and ordered from least to most broadly expressed. B. UMAP visualization of mouse lemur cell types based on the expression of hormone-related genes. Circles are cell types (unique combination of annotation name and tissue of origin) and color-coded by cell type compartment types (epithelial, neural, germ, stromal, endothelial, lymphoid, and myeloid) . Dashed lines circumscribe the cell type clusters as in panel A (also see Table 3 ), and cluster IDs and names were labeled nearby. As UMAPs qualitatively display relationships among data points, the extremely distant clusters (i.e., cluster 44. ranked from most to fewest. C. Indegree and outdegree of the network nodes (cell types), ranked from high to low values. Node degree is normalized to the total number of nodes and cluster size and measures the percentage of cell types connected from or to a network node. D. Generality score of all hormone ligands and receptors ranked from the most generally expressed to most selectively expressed. Generality score is defined as the percentage of nodes (cluster size normalized) positively expressing the gene(s) involved in ligand synthesis or receptor binding. Also see Fig. S9A for generality score defined as the number of cell type clusters positively expressing the gene(s). E. Scatterplot showing generality scores of the hormone ligands and corresponding receptors. Each circle is a unique ligand-receptor pair. F. Ranking (left) and distribution (right) of Jaccard indices of receptor pairs binding to the same ligand across all hormone receptors. Jaccard index is calculated as the number of cell types co-expressing both receptors divided by the number of cell types expressing either of the receptors. Permutation preserved node outdegree and indegree. Note that none of the permuted networks has as many 2-node feedback circuits as in the original network. C-H. Examples of two-node feedback circuits identified in the network, focusing on the endocrine cell types. Solid black arrows indicate known regulation, dashed orange arrows indicate predicted regulation with partial literature support, and dotted red arrows predicted regulation without earlier knowledge. Parentheses are used to group multiple genes that are expressed in the relevant cell type that function either together (&) or independently (/).

Concentrations of twelve hormones measured in summer (long-photoperiod) and winter (short-photoperiod) in the mouse lemur. Bars represent fold changes in hormone concentration compared to average summer levels. Error bars represent 95% confidence intervals. P-values were calculated by 2-sample t-test. Data were collected from literature and unpublished work: estradiol (Noiret et al., 2020; Perret and Aujard, 2005) , testosterone (Dal-Pan et al., 2011; Petter-Rousseaux and Picon, 1981; Schilling and Perret, 1993) , thyroxine ((Noiret et al., 2020; Petter-Rousseaux, 1984) and unpublished data), DHEA (Perret and Aujard, 2005) , cortisol ((Perret and Predine, 1984) and unpublished data), IGF1 , melatonin (Aujard et al., 2001) , and gut hormones (GIP, PPY, insulin, PYY, and GLP-1) (Dal-Pan et al., 2010; Giroud et al., 2009 ). B. Mouse lemur cell types targeted by any of these 12 seasonal hormones. Cell types are displayed in the same UMAP plot as in Fig. 2B . Circles are cell types and are color coded according to whether they express receptors for any of the 12 hormones. Dashed lines circumscribe the cell type clusters and are colored according to cell type compartments. group (red for ligands or synthases, light red for processing enzymes and/or other modulator genes, and blue for receptor genes). The dark colored dots on the intersections of the radial lines and circular lines depict mean gene expression level (by dot shade) and percent of positive cells (by dot size) for the respective gene and cell type. The outermost red or blue dots indicate positive expression of one of the hormone(s) (red) or hormone receptor(s) (blue), or both (purple) for the hormone group by the respective cell type (see Methods). The outermost grey dots indicate potentially low-quality cell types (e.g., low cell number (<10) and pancreatic cells that show notable signal contamination) whether or not the cell type expresses one the hormone ligand(s) (red outline), receptor(s) (blue outline), both (purple outline), or neither (no outline).

. Protein sequences of hormonal genes are generally more similar between human and mouse lemur than human and mouse. Each circle is a hormone-related protein. Sequences were aligned by the Needleman-Wunsch algorithm and sequence identity was calculated as the percentage of identical amino acids. Points below the red line are proteins more similar between human and mouse lemur than between human and mouse.

Circles (cell types) are color-coded by cell compartment types (epithelial, neural, germ, stromal, endothelial, lymphoid, and myeloid). Dashed lines circumscribe cell type clusters as in Fig. 2A-B and cluster IDs and names are labeled next to the cell type clusters. Cell types that were clustered by the hormone-related genes but not the non-hormonal transcriptome are labeled with asterisks (*). Note some cell types that formed discrete clusters based on hormonal gene expression intermingled with other clusters when mapped based on the non-hormonal transcriptome (e.g., cluster 51-52 and 54-58), and therefore were annotated together.

In panels A and B, cell types are displayed by the UMAP algorithm according to the cell type distances calculated with either the hormone ligands/modulators or the hormone receptors respectively. Circles (cell types) are color-coded by the cell compartment types (epithelial, neural, germ, stromal, endothelial, lymphoid, and myeloid) . Dashed lines circumscribe cell type clusters as in Fig. 2A-B and cluster IDs are annotated nearby. Cell types that were clustered using all the hormone-related genes but not in this analysis are labeled with asterisks (*). As UMAPs qualitatively display relationships among data points, the extremely distant clusters in panel A (i.e., cluster 44. mesothelial cells by hormone ligand genes) were shifted towards the center of the figure for display purposes. C. Cluster averaged silhouette values of the 68 cell type clusters by both hormone ligand and receptor genes (black), only the hormone ligand genes (red), only the hormone receptor genes (blue), or the non-hormonal transcriptome (blank). The silhouette values are a measure of how well each cell type cluster was classified (see Methods). Supplemental Tables   Table 1. Genes involved in the biosynthesis and sensing of 87 classes of hormones. Rows show hormones and columns show the hormone class name, hormone name(s), type of hormones, symbol of the genes involved in biosynthesis (Ligand/Enzyme Gene(s)) and sensing (Receptor Gene(s)), whether the hormone synthesis or maturation requires coordination of multiple tissues, years discovered, classical sites of secretion and targets, and approximate plasma concentrations in humans. Genes in the table are human genes. Rows show clusters and columns show cluster ID, name, cell types included in the cluster, and the hormones, receptors, and hormone-related genes that were positively expressed in the majority (≥50%) of the cluster cell types. Scheja, L., and Heeren, J. (2019) . The endocrine function of adipose tissues in health and cardiometabolic disease. Nat. Rev. Endocrinol. 15, 507-524. Schilling, A., and Perret, M. (1993) . Removal of the olfactory bulbs modifies the gonadal responses to photoperiod in the lesser mouse lemur (Microcebus murinus). Biol. Reprod. 49, 58-65. Schwartz, D.R., and Lazar, M.A. (2011) . Tendler, A., Bar, A., Mendelsohn-Cohen, N., Karin, O., Korem Kohanim, Y., Maimon, L., Milo, T., Raz, M., Mayo, A., Tanay, A., et al. (2021) . Hormone seasonality in medical records suggests circannual endocrine circuits. Proc. Natl. Acad. Sci. U. S. A. 118.

Terrien, J., Gaudubois, M., Champeval, D., Zaninotto, V., Roger, L., Riou, J.F., and Aujard, F. (2018) . Metabolic and genomic adaptations to winter fattening in a primate species, the grey mouse lemur (Microcebus murinus 5  10  15 20 25  30  35  40  45  50  55  60  65  70  75  7  12 13  1  28  33  37  44  47  48  51 53  57 58  62  67 68  69  71 TGFBR3  PTGER4  ADRA2B  PTH1R  PRLR  ADRA1A  FGFR4  AR  ADRA2A  OXTR  ADRA2C  THRA  SCTR  CALCRL  PAQR7  PTGFR  AVPR2  PTGER3  ADRB2  PAQR8  AVPR1A  GPRC6A  PGR  NPB  AGT  DIO1  ADM2  OXT  AVP  ACE2  FST  PTGES  IGFBP1  PTGS2  MLN  PTGS1 NPY SLC26A4 SLC40A1  ADRA1A  PRLR  AR  GCGR  PGRMC1  AVPR1A  PAQR9  ACVR1B  ITGAV  OGFR  ATP6AP2  AGT  CYP27A1  FGF21  PTGES2  HSD17B2  CFB  IGFBP1  IGFBP2  IGFBP4  KLKB1 DIO1 HAMP IGF1 

Adiponectin Receptors Form Homomers and Heteromers Exhibiting Distinct Ligand Binding and Intracellular Signaling Properties

Basic local alignment search tool

The IUPHAR/BPS Guide to PHARMACOLOGY in 2020: extending immunopharmacology content and introducing the IUPHAR/MMV Guide to MALARIA PHARMACOLOGY

Artificially accelerated aging by shortened photoperiod alters early gene expression (Fos) in the suprachiasmatic nucleus and sulfatoxymelatonin excretion in a small primate, Microcebus murinus

IGF-1: a marker of individual life-30

Estrogen receptor-alpha as a predictive biomarker in endometrioid endometrial cancer

Scale-Free and Hierarchical Structures in Complex Networks

The endothelin system and endothelin-converting enzyme in the brain: molecular and cellular studies

Molecular determinants of nephron vascular specialization in the kidney

Prolactin (PRL) and its receptor: actions, signal transduction pathways and phenotypes observed in PRL receptor knockout mice

Microcebus murinus: a useful primate model for human cerebral aging and Alzheimer's disease?

Inactivation of insulin-like factor 6 disrupts the progression of spermatogenesis at late meiotic prophase

Integrating single-cell transcriptomic data across different conditions, technologies, and species

Comprehensive single cell transcriptome lineages of a proto-vertebrate

groHMM: a computational tool for identifying unannotated and cell type-specific transcription units from global run-on sequencing data

Single-Cell RNA-Seq Technologies and Related Computational Data Analysis

Power-law distributions in empirical data

The Genetic Landscape of a Cell

Resveratrol suppresses body mass gain in a seasonal non-human primate model of obesity

Caloric restriction or resveratrol supplementation and ageing in a non-human primate: first-year outcome of the RESTRIKAL study in Microcebus murinus

Widespread seasonal gene expression reveals annual differences in human immunity and physiology

Anatomy and Physiology of the Hypothalamic-Pituitary-Gonadal (HPG) Axis. Advanced Practice in Endocrinology Nursing

Development of the endothelium: an emphasis on heterogeneity

Lipopolysaccharide-induced fever depends on prostaglandin E2 production specifically in brain endothelial cells

Mini-review: Aging of the neuroendocrine system: Insights from nonhuman primate models

The Mouse Lemur, a Genetic Model Organism for Primate Biology, Behavior, and Health

Self-perpetuating states in signal transduction: positive feedback, double-negative feedback and bistability

Adipose tissue as an endocrine organ

Insulin drives glucose-dependent insulinotropic peptide expression via glucose-dependent regulation of FoxO1 and LEF1/βcatenin

Co-expression of neuropeptide Y Y1 and Y5 receptors results in heterodimerization and altered functional properties

Mannose 6-phosphate receptors: new twists in the tale

Vasopressin Increases Urinary Acidification via V1a Receptors in Collecting Duct Intercalated Cells

A protein interaction map of Drosophila melanogaster

Gut hormones in relation to body mass and torpor pattern changes during food restriction and re-feeding in the gray mouse lemur

The Grey Mouse Lemur Uses Season-Dependent Fat or Protein Sparing Strategies to Face Chronic Food Restriction

Spermatogenesis: The Commitment to Meiosis

Vasoactive intestinal polypeptide as mediator of asthma

Circannual Rhythms: Endogenous Annual Clocks in the Organization of Seasonal Processes

The metabolic actions of glucagon revisited

Fibroblast growth factor 21 mediates specific glucagon actions

Singlecell atlas of a non-human primate reveals new pathogenic mechanisms of COVID-19

Mapping the Mouse Cell Atlas by Microwell-Seq

Paracrine signaling in islet function and survival

Stepwise loss of motilin and its specific receptor genes in rodents

Dexamethasone in Hospitalized Patients with Covid-19

Evolution of the Vertebrate Resistin Gene Family

Molecular Pathways: Clinical Applications and Future Direction of Insulinlike Growth Factor-1 Receptor Pathway Blockade

Inhibition of growth hormone signaling by the fasting-induced hormone FGF21

Endocrine FGFs: Evolution

Melanocortin crosstalk with adipose functions: ACTH directly induces insulin resistance, promotes a proinflammatory adipokine profile and stimulates UCP-1 in adipocytes

Endocrine regulation of menstruation

NCBI BLAST: a better web interface

A new model for the HPA axis explains dysregulation of stress hormones on the timescale of weeks

Handbook of Biologically Active Peptides

PGE2) promotes proliferation and invasion by enhancing SUMO-1 activity via EP4 receptor in endometrial cancer

Characterisation of the adiponectin receptors: the non-conserved N-terminal region of AdipoR2 prevents its expression at the cell-surface

A G protein-coupled receptor dimer imaging assay reveals selectively modified pharmacology of neuropeptide Y Y1/Y5 receptor heterodimers

The effects of growth hormone on adipose tissue: old observations, new mechanisms

Age-associated cerebral atrophy in mouse lemur primates

Seasonal affective disorder

The grey mouse lemur: a non-human primate model for ageing studies

Molecular evolution of NPY receptor subtypes

Resistin, a fat-derived secretory factor, promotes metastasis of MDA-MB-231 human breast cancer cells through ERM activation

The discovery of insulin revisited: lessons for the modern era

A brief history of circannual time

Progesterone activates the principal Ca2+ channel of human sperm

Functional role of progestin and the progesterone receptor in the suppression of spermatogenesis in rodents

Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets

cDNA cloning and expression of a novel adipose specific collagen-like factor

Insulin-like growth factor-1 is a negative modulator of glucagon secretion

Intrauterine Growth Retardation (IUGR) as a Novel Condition of Insulin-Like Growth Factor-1 (IGF-1)

Seasonal Variation of Overall and Cardiovascular Mortality: A Study in 19 Countries from Different Geographic Locations

Specificity and stability in topology of protein networks

UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction

Uniform Manifold Approximation and Projection (UMAP). MATLAB Central File Exchange

Seasonality of Respiratory Viral Infections

FGF21 analogs of sustained action enabled by orthogonal biosynthesis demonstrate enhanced antidiabetic pharmacology in rodents

Multifaceted roles of PGE2 in inflammation and cancer

Sex-Specific Response to Caloric Restriction After Reproductive Investment in Microcebus murinus: An Integrative Approach

Three-dimensional architecture of inner medullary vasa recta

Molecular evolutionary insights from PRLR in mammals

Human TIP49b/RUVBL2 gene: genomic structure, expression pattern, physical link to the human CGB/LHB gene cluster on chromosome 19q13

Environmental and social determinants of sexual function in the male lesser mouse lemur (Microcebus murinus)

Regulation by Photoperiod of Seasonal Changes in Body Mass and Reproductive Function in Gray Mouse Lemurs (Microcebus murinus): Differential Responses by Sex

Daily hypothermia and torpor in a tropical primate: synchronization by 24-h light-dark cycle

Aging and season affect plasma dehydroepiandrosterone sulfate (DHEA-S) 35

Effects of long-term grouping on serum cortisol levels in Microcebus murinus (Prosimii)

Age of Microcebus murinus at the onset of testicular development : Preliminary observations on photoperiodic effect

Annual variations in the plasma thyroxine level in Microcebus murinus

Annual variation in the plasma testosterone in Microcebus murinus

Promoting healthspan and lifespan with caloric restriction in primates

A draft network of ligand-receptor-mediated multicellular signalling in human

Sentence-BERT: Sentence Embeddings using Siamese BERT-Networks

THE PREPARATION, IDENTIFICATION AND ASSAY OF PROLACTIN-A HORMONE OF THE ANTERIOR PITUITARY

A proteome-scale map of the human interactome network

Asprosin, a Fasting-Induced Glucogenic Protein Hormone

Adiponectin signaling and function in insulin target tissues

Hormonal regulation of male germ cell development

Deep longitudinal multiomics profiling reveals two biological seasonal patterns in California

Physiology of proglucagon peptides: role of glucagon and GLP-1 in health and disease

Motilin: towards a new understanding of the gastrointestinal neuropharmacology and therapeutic use of motilin receptor agonists

Adipokines (Leptin, Adiponectin, Resistin) Differentially Regulate All Hormonal Cell Types in Primary Anterior Pituitary Cell Cultures from Two Primate Species

A molecular cell atlas of the human lung from single-cell RNA sequencing

Uncovering transcriptional dark matter via gene annotation independent single-cell RNA sequencing analysis

Unifying single-cell annotations based on the Cell Ontology

Hypothalamic proopiomelanocortin processing and the regulation of energy balance

Growth hormone regulation of sex-dependent liver gene expression

Deletion of prostaglandin E2 synthesizing enzymes in brain endothelial cells attenuates inflammatory fever

Evolution of the relaxin-like peptide family

The mechanism of pancreatic secretion

Cloning of adiponectin receptors that mediate antidiabetic metabolic effects

Functional and topological characterization of protein interaction networks

Neuropeptidomic Analysis Establishes a Major Role for Prohormone Convertase-2 in Neuropeptide Biosynthesis

We are grateful to Sandra Schmid, Uri Alon, Chaitan Khosla, the Krasnow lab, and the Ferrell lab for helpful discussions and comments on the manuscript. This work is supported by NIH grant R35 GM131792 to J.E. Ferrell, funding from Wall Center to M.A. Krasnow, and a fellowship from the Wu Tsai Neurosciences Institute to S. Liu. M.A. Krasnow is an investigator of the Howard Hughes Medical Institute.

 m  u  i  n  o  g  o  t  a  m  r  e  p  s e   t  y  c  o  t  a  m  r  e  p  s  y  l  r  a  e e   t  y  c  o  t  a  m  r  e  p  s  e  n  e  t  y  h  c  a  p e   t  y  c  o  t  a  m  r  e  p  s  e  n  e  t  o  l  p  i  d d   i  t  a  m  r  e  p  s  d  n  u  o  r d   i  t  a  m  r  e  p  s  g  n  i  t  a  g  n  o  l  e d   i  t  a  m  r  e  p  s  d  e  t  a  g  n  o  l  e  IGF1R  IGF2R  GALR2  PGRMC1  PAQR8  PGR  OGFR  THRA  ADRA2C  SLC40A1 UTS2R