An organism-wide atlas of hormonal signaling based on the mouse lemur single-cell transcriptome

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.


52
Hormones are chemical messengers that circulate through the bloodstream and control and 53 coordinate functions of both nearby and distant cells. Since the discovery of secretin over a 54 hundred years ago 1 , about 100 mammalian hormones have been identified 2 . Over the decades, 55 these hormones have been purified, sequenced, and synthesized, their receptors have been 56 identified, and the intracellular signaling pathways they regulate have been explored 3 . Hormones 57 regulate diverse biological processes ranging from growth and development, reproduction, 58 metabolism, and immune responses, to mood and behavior. Many hormones are also 59 therapeutically significant. For example, the discovery of insulin 4 changed type I diabetes from shown that many hormones are highly conserved across vertebrates, the hormone-producing and 72 -receiving cells are sometimes rewired for different functions 3 . 73 The grey mouse lemur (Microcebus murinus) is a small non-human primate with a lifespan of 5-74 13 years. Because of its small size, rapid reproduction, and close genetic proximity to humans, 75 the mouse lemur is an appealing model organism for studies of primate biology, behavior, aging, 76 and disease 9-13 . With the advances in single-cell RNA-sequencing (scRNAseq) techniques, we 77 have recently created a molecular cell atlas of the mouse lemur by profiling transcripts from 78 ~225,000 cells from 27 mouse lemur organs and tissues of 4 individuals 14 (https://tabula-79 microcebus.ds.czbiohub.org/). This organism-wide dataset provides an unprecedented 80 opportunity to study global hormone regulation at single cell resolution in a non-human primate. 81 In this study, we present a comprehensive analysis of primate hormone signaling based on the 82 information in the mouse lemur molecular cell atlas 14 . We systematically mapped hormone-83 producing and target cells for 87 hormone classes, and created a browsable atlas that depicts 5 GALP, and INSL4) may result from incomplete annotation of the mouse lemur genome. For 136 example, LHB encodes the beta subunit of luteinizing hormone, which is thought to be expressed 137 in the gonadotrophs of the pituitary gland across all vertebrates. The mouse lemur atlas includes 138 a small number of gonadotrophs that might be expected to express the LHB transcripts. 139 Therefore, we searched for transcripts expressed in the gonadotrophs and identified a candidate 140 gene region 20 , currently unannotated in NCBI and Ensembl, with high sequence similarity to 141 human and macaque LHB. The putative LHB region is located next to the RUVBL2 gene, as in 142 humans and mice 21 . The gene region also is selectively expressed in the expected cell type, 143 gonadotrophs ( Fig. S1-luteinizing hormone). We therefore conclude that the identified region is  Table 2). The percentages of human 152 hormone-related orthologs present in mouse lemur or mouse genomes were higher than the 153 percentage of protein-coding genes (>99.1% vs. 89.5% in the mouse lemur genome, and 98.5% 154 vs. 89.8% in the mouse genome). Most (86.3%) human hormone-related proteins were more 155 similar in sequence to their mouse lemur orthologs than their mouse homologs (Fig. S2). These 156 findings underscore the close relationship between human and mouse lemur hormone-related 157 genes.

159
Organism-wide identification of hormone-producing and target cells 160 We systematically identified the source cells and target cells for the 87 known hormone groups 161 by examining the expression of genes involved in hormone synthesis and sensing across 739 cell 162 types profiled by 10x droplet-based single cell RNA sequencing from 27 different tissues 163 included in the mouse lemur cell atlas 14 . Source cells were defined as cell types with positive 164 expression of the ligand (for protein/peptide hormones) or synthase(s) (for small molecule 165 hormones) and any necessary modulator genes. Target cells were defined as the cell types with 166 positive expression of the hormone receptor genes. We generated a browsable set of figures that 167 depict the global expression profiles of the hormone-related genes for each of the 87 hormone 168 groups as well as one set of figures for prohormone processing enzymes and one for endocrine 169 secretory proteins (Fig. S1). Prolactin (PRL). The first example is a classical endocrine hormone, prolactin. Prolactin was 180 one of the first hormones discovered 23 and has been studied in depth over the decades. Prolactin 181 is secreted by the lactotrophs of the anterior pituitary gland, the master endocrine organ. 182 Although prolactin was originally regarded as a female-specific hormone that promotes milk 183 secretion, it is now clear that prolactin plays pleiotropic roles in both sexes 2 . These roles include 184 regulation of water and electrolyte balance in the kidney, and control of skin cell proliferation 185 and hair/feather growth and shedding 24 . Prolactin also plays a major role in energy homeostasis 186 by regulating cell growth and metabolism in the liver, adipose tissues, and pancreas in both 187 sexes 25 , and may be involved in the immune response 26 . These sex-non-specific functions are 188 conserved, not just among milk-producing mammals, but across all vertebrates. 189 As expected, prolactin transcripts were selectively expressed in the lactotrophs of the mouse 190 lemur pituitary gland, and were not detected in other tissues (Fig. 1, Fig. S1-prolactin). The 191 prolactin receptor (PRLR) was found in both male and female lemurs in the expected cell types, 192 including kidney nephron tubule cells, skin epithelial cells, pancreatic ductal and acinar cells, 193 hepatocytes, intestinal epithelial cells, adipocytes, reticular cells, plasma cells, choroid plexus 194 cells, and the neighboring TSH-secreting thyrotrophs and FSH/LH-secreting gonadotrophs in the 195 pituitary gland (Fig. 1, Fig. S1-prolactin). Thus, single-cell RNA sequencing confirms that 196 prolactin signaling is similar in mouse lemurs and other animals. 198 in body weight and fat reserves [27][28][29] . Adipocytes were originally regarded only as nutrient storage 199 cells, but it is now clear that they also serve an endocrine role 30 . Adiponectin and leptin are two 200 important metabolic hormones secreted by adipocytes 31,32 . Surprisingly, leptin transcripts were 201 almost undetectable in the adipocytes captured in the atlas (see detailed discussion in the 202 accompanying Tabula Microcebus manuscript 2) ( Fig. S1-leptin). 203 However, adiponectin was abundantly and specifically expressed in mouse lemur adipocytes 204 captured from multiple tissues, including four adipose depots (subcutaneous, mesenteric, 205 interscapular brown, and perigonadal), as well as kidney, limb muscle, pancreas, and bone ( Fig.   206 1, Fig. S1-adiponectin). Both UCP1-low white-like and UCP1-high brown-like adipocytes 207 highly expressed ADIPOQ. 208 Previous research has characterized tissue-specific functions of adiponectin in the liver, 209 pancreas, skeletal muscle, endothelial cells, macrophages, and adipocytes 33 . While we indeed 210 observed expression of the adiponectin receptor ADIPOR1 in these tissues in the mouse lemur, 211 we also found ADIPOR1 to be abundantly expressed in almost all lemur cell types (Fig. 1). The 212 ubiquitous expression of this receptor suggests that adiponectin exerts an organism-wide effect 213 through ADIPOR1. 214 In comparison, ADIPOR2, the other adiponectin receptor, was expressed more selectively, most 215 notably in the male germ cells and adipocytes. Interestingly, the mouse lemur expression pattern 216 of ADIPOR2 differs from that of the mouse in which ADIPOR2 is reported to be most highly 217 expressed in the liver and almost undetectable in the testis 34 , suggesting differential ADIPOR2 218 effects in the two species. 219 Motilin (MLN). Comparing expression patterns of the mouse lemur to that of the mouse and 220 human, we identified interesting cases of hormone network rewiring. A telling example is 221 7 motilin, canonically considered to be a gastrointestinal hormone secreted by enteroendocrine 222 cells to regulate gastrointestinal smooth muscle movement 15 . Mice lack both the motilin and the 223 motilin receptor genes (Table 2). Indeed, genome comparative studies reported that motilin 224 signaling is conserved in vertebrates but is lost in rodents 35 . The lack of motilin signaling in mice 225 and rats has hampered progress on further characterizing motilin function as well as the 226 development of motilin agonists for the treatment of gastroparesis 36 . In the mouse lemur, we 227 observed motilin expression in enteroendocrine cells in the small intestine but not colon ( Fig.   228 1A), which is consistent with the canonical view. However, there was also a broad range of other 229 motilin-expressing cell types, including kidney nephron tubule cells, chemosensory tuft cells 230 from both the airway (lung and trachea) and small intestine, certain brain neurons, adipocytes 231 (higher in UCP1-high brown-like adipocytes than in white-like adipocytes), mesothelial cells, 232 and the majority of the CD4+CD8-CCL5-helper T cells (Fig. 1, Fig. S1-motilin). The receptor 233 was found to be most highly expressed in pancreatic delta cells and several epithelial cell types 234 of the skin (e.g., sweat gland cells). Together, this expression pattern suggests that motilin is 235 involved in a diverse array of mouse lemur physiological functions in addition to the regulation 236 of gut motility.

237
Resistin (RETN). We also noted interesting species differences in resistin signaling. Resistin is 238 a phylogenetically recent hormone, present only in mammals 2 . It is expressed mainly by 239 adipocytes in mice but by macrophages, monocytes, and neutrophils in humans 37,38 . In the mouse 240 lemur, we found that resistin is highly expressed in macrophages, monocytes, and neutrophils, 241 and nearly absent from adipocytes ( Fig. 1, Fig. S1-resistin), as it is in humans. Interestingly, 242 resistin is also highly expressed in many of the vascular endothelial cell types and in metastatic 243 uterine tumor cells, suggesting additional sources for this hormone. 244 Prohormone processing enzymes. Prohormone processing enzymes are required for the 245 maturation of many peptide hormones, including insulin, adrenocorticotropic hormone (ACTH), 246 gastric inhibitory polypeptide (GIP), and somatostatin 15 . The required processing enzymes have 247 been identified for some but not all hormones 39,40 . In this study, we considered both the ligand 248 and the known or suggested processing enzymes when determining if a cell type was a hormone 249 source. Selective expression of the processing enzymes can help distinguish hormones that share 250 the same precursor gene. For example, proglucagon (GCG) is the precursor for both glucagon 251 and glucagon-like peptide-1 (GLP-1). GCG was expressed in both pancreatic α cells and 252 intestinal enteroendocrine cells, but the two cell types expressed different processing enzymes: α 253 cells expressed only PCSK2 which is required for glucagon synthesis, and enteroendocrine cells 254 expressed only PCSK1 which is required for GLP-1 synthesis (Fig. S1-glucagon). Thus, as is the 255 case in mice and humans 41 , α cells are sources of glucagon and enteroendocrine cells are sources 256 of GLP-1 in the mouse lemur. 257 Interestingly, we noticed that some cell types in the mouse lemur transcribed prohormone genes 258 without expressing the respective processing enzymes. For example, proopiomelanocortin 259 (POMC) is the precursor of multiple hormones including ACTH, melanocyte-stimulating 260 hormone (MSH), β-endorphin, and lipotropin, and the first step in POMC processing is carried 261 out by PCSK1 42 . As expected, ACTH-secreting corticotrophs expressed both POMC and the 262 necessary processing enzymes, but mouse lemur male germ cells expressed abundant POMC 263 without PCSK1 (Fig. S1-proopiomelanocortin-derived hormones). Similarly, male germ cells 264 expressed GIP, which also requires PCSK1 processing. We also found renin (REN) expression in  These results suggest several possible hypotheses: different peptidases are involved, the   269   prohormone transcripts are not translated, the prohormone proteins may function without   270 proteolysis, or the processing enzymes may be expressed below detection levels in the animals 271 that were sampled. 272 In summary, we have depicted the organism-wide, single-cell expression patterns for hormones, 273 synthases, hormone modulators, and receptors from 87 hormone groups, and have systematically 274 identified their sources and targets (Fig. S1). This global analysis detected canonical cell types Hormone-related gene expression faithfully classifies cell type 281 Next, we clustered mouse lemur cell types based on the transcriptional profiles of only the 282 hormone-related genes ( Fig. 2A-B). We organized the cell types into 75 clusters and summarized 283 the hormones and hormone receptors that were expressed in the majority of the cell types for 284 each cluster in Table 3. We noted that cells of the same type almost always clustered closely 285 irrespective of their tissue of origin, and different cell types were well-separated in the cluster 286 dendrogram ( Fig. 2A). For example, all lymphatic endothelial cells clustered closely and showed 287 characteristic IGF1 expression despite their coming from nine different tissues and four 288 individuals (cluster 53, Fig. 3B, Fig. S1-insulin-like growth factor, Table 3). Schwann cells 289 from nine tissues clustered into myelinating and non-myelinating subtypes (cluster 40 and 37). 290 Adipocytes sampled from four adipose depots and six other tissues formed one cluster (cluster 291 13, Table 3). Skeletal muscle stem cells from the limb muscle and tongue clustered together 292 (cluster 17, Table 3). Mesothelial cells sampled from eight visceral organs formed a cluster 293 distant from other stromal cells and showed distinctive expression of NPR1, CFB, PTGIS, 294 ITGB3, and PNMT (cluster 44, Fig. 2A-B, Fig. S4A, Table 3). Although the analysis only 295 examined the transcripts of 329 hormone-related genes, ~1% of all detected genes, these results 296 suggest that hormone-related gene expression is indicative of the overall transcriptome among 297 different cell types. 298 To further test this hypothesis, three comparative analyses were performed. First, we projected 299 the high-dimensional transcriptional profiles of the mouse lemur cell types into two-dimensional 300 uniform manifold approximation and projection (UMAP) representations, which illustrate cell 301 type similarity in the space of the selected genes. We compared the spatial relationship of cell 302 types in two UMAPs: one constructed based on only the 329 hormone-related genes (Fig. 2B), 303 the other based on all genes excluding the hormone genes (>30,000 genes) (Fig. S3). In each 304 case, the data points represent cell types rather than individual cells (Fig. 2B, Fig. S3, Fig. S4). 305 The two UMAPs showed remarkable similarity and cell types were similarly clustered with only 306 minor differences (see Fig. S4C for comparison). One such difference was the chemosensory tuft 307 cells that reside in the epithelium of both the airway (lung, trachea) and small intestine. In the 308 9 space of hormone-related genes, tuft cells from all three tissues formed a cluster (cluster 21) with 309 distinctive expression of renin (REN), cyclooxygenase 1 (PTGS1), arachidonate 5-lipoxygenase- 310 activating protein (ALOX5AP, required for leukotriene production), and motilin (MLN) (Fig.   311 2B, Fig. S1, Table 3). However, in the space of all genes, the lung, trachea, and intestinal tuft 312 cells were separated from each other (labeled as 21*A, 21*B, and 21*C in Fig. S3). Another 313 notable exception was skin melanocytes (cluster 38), which clustered closely with 314 oligodendrocytes (cluster 39) when using only the hormone-related genes ( Fig. 2B) but clustered 315 with skin epithelial cell types (cluster 36) in the space of all genes (Fig. S3). These differences 316 notwithstanding, the hormone-gene and transcriptome-based cell type UMAPs were qualitatively 317 similar and shared an overall clustering pattern. 318 We then compared the transcriptional distances between pairs of cell types for either just the 319 hormone-related genes (hormone-based distances) or the non-hormonal transcriptome 320 (transcriptome-based distances) (Fig. 2C, Fig. S5A-B) Ontology graph and the text-based distances from the text descriptions of the Cell Ontology 326 terms 44 . These two distances were also significantly correlated, with r=0.45 (p < 1E-15) ( Fig.   327 2D, Fig. S5A and C). Taken together, these analyses support the notion that the small number of 328 hormone-related genes faithfully classifies most of the cell types of the mouse lemur. 329 We then examined if the expression of hormone ligands or receptors alone could distinguish cell 330 types (Fig. S4). Interestingly, for most cell type clusters, hormone ligands or receptors alone  subtypes (hepatocytes), and primary vs. metastatic tumor cells (Fig. 3).

10
Spermatogenesis. Spermatogenesis is a tightly regulated process that happens asynchronously 351 and continuously in the male testis 45 . Therefore, sampling the testis at one fixed time point can 352 capture snapshots of the entire process. The atlas includes seven sperm and sperm progenitor cell 353 types that are at different spermatogenesis stages as indicated by the expression of known stage 354 markers 14 . Interestingly, these cells exhibited spermatogenesis-stage-dependent hormone 355 signaling (Fig. 3A). The seven cell types fall into three major stages: the undifferentiated 356 spermatogonia, the spermatocytes (which undergo meiosis), and the haploid spermatids that 357 elongate and form mature sperms. Comparing hormone-related gene expression among these 358 three major stages (cluster 1, 2, 8), over 50 hormone ligands, enzymes, and receptors were 359 differentially expressed (Fig. 3A), suggesting an extensive shift in hormone signaling. Previous 360 studies on spermatogenesis have mainly focused on testosterone and on the supporting cells 361 rather than sperms and sperm progenitors 46,47 . The present data suggest that sperms and sperm 362 progenitor cells are not only differentially regulated by a variety of hormones but may also 363 secrete a changing array of signaling molecules as they mature. Some of these molecules have 364 been recently suggested to be involved in sperm functions. For example, progesterone was 365 known to activate a sperm-specific calcium channel prior fertilization to induce Ca 2+ influx, but 366 the underlying progesterone receptor(s) remained unidentified 48 . Two non-classical progesterone 367 membrane receptors (PAQR5/9) were highly expressed in the lemur spermatids and therefore 368 may mediate the process (Fig. 3A). Interestingly, two additional membrane progesterone 369 receptors (PAQR8, PGRMC1) were specifically expressed in the undifferentiated spermatogonia 370 (Fig. 3A), suggesting that they may contribute to regulation of early spermatogenesis by 371 progesterone 49 . Insulin-like factor 6 (INSL6), a member of the secretory relaxin family, was 372 selectively expressed in the lemur spermatocytes (Fig. 3A). Supporting a function in these cells,  Fig. 2A-B, Fig. S3), the lung 381 and brain vascular endothelial cells exhibited unique hormone signaling signatures (Fig. 3B). In contrast, brain vascular endothelium was found to express higher levels of a few hormone 399 ligands/synthases. Endothelin-3 (ET-3), a member of the vasoconstricting endothelin family, was 400 selectively expressed in the brain vascular endothelium (Fig. 3B). ET-3 peptide has been 401 detected at various brain locations in rodents, but it remained unsettled what cells secreted the 402 peptide and what physiological role it played 55 . Our results suggest that ET-3 is synthesized by 403 the endothelial cells in the brain, which could target the nearby EDNRB-expressing SLC7A10+ 404 astrocytes and the EDNRA-expressing vascular smooth muscle cells and pericytes of the brain 405 ( Fig. S1-endothelin). Additionally, the brain endothelium also expressed high levels of PTGES 406 and PTGS2/COX-2 (Fig. 3B), which are the synthases for prostaglandin E2 (PGE2). PGE2 has 407 been suggested to regulate inflammation-induced fever in mice 56,57 . Interestingly, brain 408 endothelial cells also selectively lacked expression of IGFBP4, which was abundantly expressed 409 in the endothelium of all other tissues (Fig. 3B). 410 The expression of the detected hormone signature genes (RAMP1, CALCRL, and VIPR1 in 411 lung, and EDN3, PTGES, and PTGS2 in the brain) was higher in endothelial cells than other cell 412 types from the same organ (Fig. S1). Thus, it is unlikely that these organ-specific endothelial 413 expression patterns resulted from cross-contamination (signal spreading). Altogether, lung and 414 brain endothelial cells exhibited organ-specific hormone signaling which likely regulates organ-415 specific physiological functions.

421
Nephrons are single-cell layered, one-directional tubules that serve as the structural and 422 functional units of the kidney. We identified a large number of hormone-related genes that are 423 specific to different sections of the nephron epithelium (Fig. 3C). Some of these hormones target 424 specific nephron sections to regulate urine formation. For example, vasopressin is well known 425 for its antidiuretic effect, promoting water reabsorption in the collecting duct principal cells 2 . 426 Indeed, the vasopressin receptor AVPR2 was expressed exclusively in the collecting duct 427 principal cells. Interestingly, the collecting duct intercalated cells expressed another vasopressin 428 receptor, AVPR1A (Fig. 3C), which has been recently reported to mediate the regulation of acid 429 homeostasis by vasopressin 58 . Although the kidney epithelium is not commonly considered to be 430 an endocrine tissue, many hormone ligands/synthases were also expressed in different nephron 431 sections, including angiotensin (AGT), adrenomedullin 2 (ADM2), motilin (MLN), and 432 neuropeptide Y (NPY). This suggests that nephron epithelium may play multifaceted roles in 433 hormone production in addition to urine formation. 434 We also identified interesting differences in hormone signaling in the vascular endothelium of 435 different regions of the mouse lemur kidney. In particular, we identified a gradual shift in 436 hormone signaling along the arteriole-to-venule-axis of the vasa recta endothelium (Fig. 3D). Interestingly, the metastatic tumor cells also acquired new abilities in hormone production. For 484 example, elevated expression of prostaglandin E2 synthase (PTGES2) was present in metastatic 485 cells (Fig. 3F). PGE2 has been suggested to be a tumorigenic factor in many cancers 64 . A recent 486 study identified that the PGE2 synthase PTGES2 was elevated in human endometrial cancer 487 tissues and found that PGE2 promotes cell proliferation and invasion in in vitro cell assays 65 . 488 Resistin was also found to be highly expressed in the metastatic cells (Fig. 3F) The hormone network was remarkably densely connected; 51% of all possible directed edges 509 were found to be present. In fact, in the graphical depiction of the network (Fig. 4A), individual 510 edges are hardly distinguishable because the edges are so dense. The density of the hormone 511 network is comparable to that of the cytokine network 67 , but much higher than the densities of 512 other characterized biological networks. For example, the network densities of the yeast, fly, and 513 human protein-protein interaction networks are all estimated to be ~0.1% 68-70 , and the yeast 514 genetic interaction network density is ~3% 71 . The dense connection between cell types by 515 hormonal regulation is also highly robust: 67% of all existing edges are connected by more than 516 one ligand-receptor pair, and the network maintains high density connection when randomly 517 removing 10% of the ligand-receptor pairs (Fig. S7A). 518 The density of the hormone network is largely due to the fact that multiple hormone ligands and 519 receptors are expressed in most cell types (Fig. 4B). On average, cell types expressed 17.0 ± 7.9 520 (mean ± S.D.) receptors or receptor complexes. In comparison, fewer hormone ligands were 521 14 expressed (2.8 ± 2.5, mean ± S.D.). Only a few cell types expressed 10 or more hormone ligands 522 (e.g., enteroendocrine cells, pancreatic endocrine cells, certain neurons, and certain fibroblasts), 523 and 38.3% of cell types expressed a single ligand or no ligand. This indicates that cells in general 524 respond to more hormones than they produce. 525 The five neuroendocrine cell types of the anterior pituitary-corticotrophs, thyrotrophs, 526 gonadotrophs, lactotrophs, and somatotrophs-play a vital role in endocrinology. These cells are 527 components of regulatory cascades, and classically they are regarded as receiving regulatory 528 signals from one or a few hormones produced in the hypothalamus, and in turn producing one or 529 a few hormones to regulate a peripheral endocrine gland. Yet even these cells were found to 530 express numerous receptors, ranging from 17 (corticotrophs) to 29 (gonadotrophs and 531 somatotrophs) ( Table 4). These cells appeared to be regulated by growth factors, neuropeptides, 532 peripheral and classical hypothalamic hormones, as well as neurotransmitters (e.g., dopamine) 533 that were not on our hormone list. Likewise, they were found to express genes for several 534 hormones in addition to the classical hormones (Table 4) Table   553 1). 554 Some ligands and receptors were expressed in highly specific cell types while others were more 555 broadly expressed. We therefore quantified for each ligand or receptor a generality score to 556 evaluate the range of its expression across different cell types in the atlas ( Fig. 4D-E, Fig. S9, 557 Table 5). Strikingly, a small number of hormone receptors were globally expressed. For (69.5%), the opioid growth factor receptor OGFR (65.9%), and the thyroid hormone receptor 563 THRA in complex with RXRB (54.3%). While it has been recognized that some of these 564 hormones exert organism-wide effects, our results suggest that they target not only a broad range 565 of tissues but likely regulate a large portion or perhaps majority of the cells of the body. 566 In comparison, most hormone ligands and receptors were expressed in specific cell types; ~85% 567 of the hormone ligands and ~60% of the hormone receptors were detected in less than 5% of the 568 cell types. Certain hormones like follicle stimulating hormone (FSHB and CGA) were only 569 detected in a single cell type. Only a few hormone ligands were relatively broadly expressed, 570 including adrenomedullin (ADM, 18.6%), acylation stimulating protein (C3, 18.1%), IGF2 571 (17.9%), and resistin (RETN, 15.4%). The generality score followed a power-law distribution 572 toward the tail where a small number of globally expressed receptors dominated the distribution 573 ( Fig. S9B-C). This resulted in a bimodal distribution of the network node outdegree (abundance 574 of outgoing edges/downstream regulation) ( Fig. 4A and C, Fig. S8). Cells expressing hormone 575 ligand(s) to the high-generality receptors (e.g., adipocytes which secrete adiponectin) sent 576 globally connected outgoing edges, and other cells that only expressed ligand(s) to specific 577 receptors sent a much smaller number of outgoing edges. In contrast, the nodes' indegree 578 (abundance of incoming edges/upstream regulation) was more homogeneously distributed 579 because of the higher specificity of hormone ligand expression and large number of hormone 580 receptors expressed in the cells (Fig. 4A-C, Fig. S8).

582
Receptors for the same hormone tend to be expressed in mutually exclusive cell types 583 Many hormones have more than one receptor, which may mediate different functions ( Table 1). 584 We therefore examined if receptors for the same ligand tend to be co-expressed or if their 585 expression is mutually exclusive across different cell types. To quantitatively compare the target 586 cell type profiles, we calculated a Jaccard index (JI, intersection over union) for each qualifying 587 receptor pair, namely the ratio of cell types that co-expressed both receptors to the cell types that 588 expressed either of the receptors. This analysis revealed that the majority of the receptor pairs 589 had extremely low JIs: 25% (41/161) of the receptor pairs had a JI equaling zero, meaning no 590 overlap in their target cell types, and 57% (92/161) of the receptor pairs had a JI of less than 5% 591 (Fig. 4F). A telling case is the adrenoceptors; although epinephrine (adrenaline) targets many 592 cell types, different adrenoceptor subtypes were selectively employed in different cell types ( Fig.   593 4F, Fig. S1-epinephrine/norepinephrine). For example, adrenoceptors ADRA1D and ADRB3 594 shared no target cell types; ADRA1D was expressed broadly in vascular smooth muscle cells and 595 pericytes of multiple organs, certain glial cells, and leptomeningeal cells, whereas ADRB3 596 mainly targeted adipocytes. Vasopressin receptors AVPR1A and AVPR2 also had less than 3% 597 overlap in their target cell types; AVPR2 was predominantly expressed in the kidney collecting 598 duct principal cells, while AVPR1A was highly expressed in the kidney collecting duct 599 intercalated cells and broadly expressed in many stromal cell types (Fig. S1-vasopressin). Hormone regulation often involves feedback control. We therefore searched the mouse lemur 618 hormone cell communication network to systematically identify potential feedback regulation. 619 For simplicity, we focused on two-node feedback circuits, namely cell type pairs that are 620 connected by edges in both directions (Fig. 5A). This systematic search revealed a vast number 621 of feedback circuits; 29% of all node pairs formed feedback circuits, showing that feedback 622 control is widespread in the hormone network ( Table 6). We then tested if feedback circuits are 623 enriched in the hormone network, i.e. more prevalent than would be expected by chance. We 624 applied degree-preserving random rewiring 80 to randomly permute edge connections without 625 changing node degree (i.e., network edge number and distribution). Fig. 5B shows the 626 distribution of the feedback circuit counts from such randomly permuted networks (n = 1000). 627 The number of feedback circuits in the permuted networks was consistently lower than in the 628 original hormone network. However, the difference was small (on average, only ~0.2%). 629 Additionally, we randomly removed 10% of ligand-receptor pairs to test how perturbations to 630 ligand-receptor binding affect feedback distribution of the network. This analysis revealed that 631 the network feedback circuits remained abundant despite edge removal as cell connections 632 usually involved multiple ligand-receptor pairs (Fig. S7B, Fig. 4B). Taken together, this analysis 633 suggests that feedback circuits are abundantly and robustly distributed as a result of the dense 634 connection of the hormone network. Evolution may have further enriched feedback circuits in 635 hormone regulation as suggested by the random rewiring analysis. 636 We next reviewed selected feedback circuits between endocrine cell types ( Fig. 5C-L hormones, which feed back to gonadotrophs through androgen and estrogen receptors (AR, 647 ESR1). Additional feedback regulation from granulosa cells to gonadotrophs was also mediated 648 through anti-Müllerian hormone (AMH) and inhibin/activin signaling. We also detected the well-649 established feedback loop between somatotrophs in the anterior pituitary and hepatocytes in the 650 liver, which was mediated through growth hormone (GH), IGF1 2,82 , and potentially hepcidin and 651 FGF21 signaling 83 (Fig. 5D). Additionally, in the pancreas, we detected the known feedback 652 circuits between pancreatic α and β cells through glucagon and insulin signaling and between α 653 and δ cells through glucagon and somatostatin signaling (Fig. 5E-F) 84 . The detection of these 654 various known feedback loops supports the validity of the comprehensive scRNAseq approach. 655 The analysis also revealed a long list of potential endocrine feedback circuits that have not been is specific to primates. 680 We also detected feedback circuits of intestinal enteroendocrine cells with pancreatic β cells and 681 adipocytes ( Fig. 5K-L). It is known that enteroendocrine cells secrete GIP, which induces reproductive organs for mating. The global nature of these phenotypic changes suggests that they 702 may be hormonally mediated. To date, only 12 circulating hormones have been measured in 703 mouse lemurs, but all showed dramatic seasonal changes (Fig. 6A, Table 7). This suggests that 704 hormones are likely the key messengers coordinating and driving various seasonal changes 705 among different organs. The hormone atlas suggests that these seasonally changing hormones 706 target a broad array of cells and tissues; 95% of the cell types possessed receptors for at least one 707 of these 12 hormones. This finding suggests that seasonal rhythms engage not only tissues that 708 demonstrate obvious morphological changes (e.g., adipose tissues and gonads), but involve 709 global shifts in cellular metabolism and physiology (Fig. 6B).

711
Hormonal signaling is vital to multicellular life, allowing cells to communicate with both 712 surrounding cells and remote cells in other tissues. In this study, we systematically mapped 713 hormone-producing and -target cells across the 739 cell types from 27 different tissues from the 714 mouse lemur molecular cell atlas. We generated a browsable atlas that depicts gene expression of 715 ligands, modulators, and receptors for 87 hormone classes (Fig. S1, Fig. 1 less than 1% of the transcriptome (Fig. 2, Fig. S3, Fig. S5). Ligands alone and receptors alone 724 also clustered cell types well, with certain clusters mostly defined by their ligands and others by 725 their receptors (Fig. S4). It is remarkable, perhaps, how well a cell's intercellular communication 726 repertoire encapsulates the cell's identity. 727 There were, however, some exceptions to this general rule. For example, sperms and sperm 728 progenitor cells exhibited dramatic changes in hormone gene expression during spermatogenesis 729 19 despite their close developmental connection (Fig. 3A). In addition, vascular endothelial cells 730 from lung and brain exhibited hormone signatures that were distinct from those of other vascular 731 endothelial cells (Fig. 3B). Hormone signaling varied in an orderly and graded fashion along the 732 epithelium and vasa recta endothelium of the nephron (Fig. 3C and D). Hepatocytes were found 733 to fall into two molecular subtypes with differential hormonal gene expression (Fig. 3E, Fig.   734 S6), suggesting different physiological roles. And finally, metastatic endometrial tumor cells 735 exhibited different hormonal signaling from that seen in their likely cells of origin in the uterus 736 (Fig. 3F), suggesting possible rewiring of hormone signaling during tumor development and 737 metastasis. 738 The cell-type specific gene expression data allowed us to construct an organism-wide cell 739 communication network for mouse lemur hormone signaling, by connecting the hormone-740 producing cells to the corresponding target cells (Fig. 4). This hormone network was 741 exceptionally dense and robust to deletion of individual edges (Fig. 4A, Fig. S7A), which mainly 742 resulted from two network features. First, most cells in the network expressed multiple hormone 743 receptors and/or ligands, though in general cells were regulated by more hormones than they 744 produced (Fig. 4B). Second, although most hormone ligands and receptors were expressed in 745 specific cell types, a small fraction of hormone receptors were nearly ubiquitously expressed 746 (e.g., adiponectin receptors, a membrane progesterone receptor, the IGF1 receptor, the thyroid 747 hormone receptor, and the cortisol receptor), resulting in almost global connection of the relevant 748 hormone-producing cells to the other cells of the body (Fig. 4C-E, Fig. S8, Fig. S9). 749 Interestingly, despite the density of the hormone signaling network, subtypes of hormone 750 receptors that bind to the same ligand were generally expressed in a mutually exclusive fashion, 751 supporting the notion that different hormone receptor subtypes usually mediate different 752 physiological functions through separate target cells (Fig. 4F) Feedback regulation was widespread in hormonal signaling, and it was robust to network 759 perturbations (Fig. 5B, Fig. S7B). The analysis identified well-known feedback circuits among 760 several endocrine cell pairs and detected additional, previously unreported, hormones involved in 761 feedback control (Fig. 5C-F). We also identified many novel feedback circuits (Fig. 5G- (Fig. 6A). The hormone atlas shows that almost all (95%) of the mouse lemur cell types are 771 regulated by seasonally varying hormones (Fig. 6B). This suggests that most, if not all, mouse were first normalized by total transcript counts, multiplied by 10,000, averaged across all cells in 849 the cell type, then natural log transformed (i.e., log(1+AverageUMIcountPer10KTranscripts); 850 note that 1 is added to the average so that genes with zero transcripts will yield zero rather than 851 negative infinity). To determine positive expression at the cell type level, we first applied a low-852 level absolute threshold on average expression level (0.037) and percent of nonzero cells (5%) 853 below either of which we considered the expression to be negative. Second, we calculated a 854 "maximal" expression level of the gene by searching the high-expression levels across the atlas. 855 Specifically, we considered the 99th percentile of the cells that passed the first thresholding to be 856 the maximum to increase the robustness of the analysis. Lastly, we applied relative thresholding 857 of the "maximal" level (for ligand genes at 20% and other gene types at 5%). Cell types that 858 22 passed all thresholding criteria were considered to have functional expression of the gene. We 859 used a higher relative threshold for the ligand genes because 1) ligands are secreted into the 860 circulation so a highly expressing source outshadows a low expressing source, whereas non-861 ligand products of other genes (e.g., receptors, synthases, processing enzymes) function inside 862 individual cells; 2) we noted that ligand genes were usually expressed at higher levels and were 863 more easily affected by cross-contamination (signal spreading) inside a tissue. For hormone 864 ligands and receptors that included multiple genes, we applied AND or OR logic gates as 865 indicated in Table 1 to determine if a cell type expressed the ligand or receptor. For an AND 866 logic gate, a functional cell type should express all its components above the functional level; 867 while for an OR logic gate, expression of any of its components would suffice.  Table 3 with the list of hormones, hormone receptors, and hormone-related genes that 894 were expressed in the majority (≥ 50%) of the cell types for each cluster. 895 In addition to visualizing the cell type clustering via dendrogram, we also applied uniform 896 manifold approximation and projection (UMAP) 118 to embed the high dimensional expression 897 data to 2D using a Matlab UMAP package 119 . Cell type pairwise cosine distances were used as 898 input when generating the UMAPs. 899 For comparative analysis, we additionally calculated distances between cell types by either the tissues and would be categorized as different "cell types" by our definition. We did not merge 978 circulating cell types from different tissues because it is challenging for some cell types like T 979 cells to determine whether they are circulating or potentially tissue specific or tissue-resident. 980 We also tried simply normalizing the node degree by the total number of nodes; similar 981 conclusions were reached using either method. Calculating target cell Jaccard indices for different receptors that bind to the same ligand: 1001 We calculated Jaccard indices for pairs of receptors that 1) can bind to the same ligand as 1002 described in Table 1   additional data in Table 7), IGF1 130 , melatonin 131 , and gut hormones (GIP, PPY 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 NPR2  FGFR1  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