A molecular map of murine lymph node blood vascular endothelium at single cell resolution

Blood vascular endothelial cells (BECs) control the immune response by regulating blood flow and immune cell recruitment in lymphoid tissues. However, the diversity of BEC and their origins during immune angiogenesis remain unclear. Here we profile transcriptomes of BEC from peripheral lymph nodes and map phenotypes to the vasculature. We identify multiple subsets, including a medullary venous population whose gene signature predicts a selective role in myeloid cell (vs lymphocyte) recruitment to the medulla, confirmed by videomicroscopy. We define five capillary subsets, including a capillary resident precursor (CRP) that displays stem cell and migratory gene signatures, and contributes to homeostatic BEC turnover and to neogenesis of high endothelium after immunization. Cell alignments show retention of developmental programs along trajectories from CRP to mature venous and arterial populations. Our single cell atlas provides a molecular roadmap of the lymph node blood vasculature and defines subset specialization for leukocyte recruitment and vascular homeostasis.

T he vascular endothelium lining blood vessels regulates exchange of oxygen and metabolites between the blood vascular compartment and tissues. In lymph nodes (LN), additionally, the endothelium plays essential roles in controlling immune cell access. The organization of the vasculature in LN is well characterized: Arteries entering at the hilus lead to capillary arcades in the LN cortex; capillaries link to 'high endothelial venules' (HEV), the post capillary venules that recruit lymphocytes from the blood 1 ; and the vasculature exits the lymph node at the hilus. Upon immunization, lymph nodes can increase in volume 10-fold or more within days, and the vascular endothelium expands roughly proportionally from local precursors 1,2 . Vessel expansion involves extensive proliferation of both capillary and high endothelial cells (HEC) 3 , without contribution from blood borne progenitors 4 ; and an increase in vessel numbers through intussusceptive (splitting) angiogenesis 5 . However, the nature and extent of endothelial cell diversity within LN remains incompletely understood.
Single-cell RNA profiling is a transformative technology for the identification of cell diversity and elucidation of developmental and physical relationships. Here we provide a survey of blood vessel endothelial cells (BEC) from mouse peripheral lymph nodes (PLN). We confirm known features of high endothelium 6,7 , identify endothelial subsets, describe unexpected diversity among capillary cells and demonstrate a distinctive role of medullary veins in selective myeloid cell recruitment. We define gene signatures and transcription regulatory factors for these subsets, and map key subsets to the vasculature. We also identify a primed capillary resident regenerative population (CRP) that displays the angiogenic endothelial marker Apln, is enriched in cells undergoing cell division, and possesses stem cell and migratory gene signatures. Genetic lineage tracing suggests that CRP function as progenitors, contributing to neogenesis of the blood vascular endothelium including high endothelium in immune angiogenesis.

Results
Single-cell profiling of lymph node blood endothelial cells. We performed single-cell RNA sequencing (scRNAseq, 10× Chromium) of sorted BEC from PLN of adult mice (Fig. 1a). Four cohorts were analyzed including a group of male and female Balb/c mice (PLN1) processed together and resolved in silico (PLN1_m and PLN1_f), an additional group of female Balb/c mice (PLN2), and one of female mice of a mixed background (PLN3) (Supplementary Fig. 1a). Unsupervised analyses (Methods) defined 8 BEC subsets with distinct gene expression: arterial EC (Art); 2 venous EC subsets, high endothelial cells (HEC) and non-HEC veins (Vn); and five capillary phenotype EC subsets, including a transitional phenotype capillary EC (TrEC) and primed EC comprising a capillary resident regenerative population (CRP) (Fig. 1b-d). Annotations of the arterial and venous populations were guided by known gene markers 6,8,9 . Nearest neighbor alignments, which can model developmental relationships or spatial alignments of cells, were visualized by trajectory inference using tSpace 10 , which revealed a continuum of cell phenotypes with branching of arterial and venous subsets from capillary EC populations (Fig. 1e). Clusters extending to the termini of the arterial, CRP, and HEV branches were further separated to distinguish cells that were most distinct (distant along trajectories) from the bulk of EC (darker shading, Fig. 1e): these cells are enriched for genes associated with mature arterial (e.g., Gkn3, Bmx) or HEV (Chst4, Glycam1) differentiation, or in the case of CRP, for genes associated with endothelial specification during early development or with stem cells (e.g., Ets2, Cxcr4, Nes) (Fig. 1f). Clusters and trajectory alignments were shared by male and female mice and by the independently processed samples ( Supplementary Fig. 1b-d). Gene expression signatures were generated for each subset using the combined scRNAseq datasets (Methods; Supplementary Data 1). Correlation in mean gene expression profiles of the identified subsets across the cohorts is shown in Supplementary Fig. 1e. Expression of marker genes depicted in Fig. 1f showed overall consistency across all biological replicates ( Supplementary Fig. 2).
Characterization of the arterial and venous subsets. The arterial cluster among the profiled BEC was identified by expression of Gja5 and Gja4 (encoding connexin 37 and 40, respectively) and Bmx (Fig. 1f). Consistent with prior reports describing Gkn3 as a marker for mature arteries 8 , it is selectively expressed in mature Art, and absent in pre-Art, which lay closer to the capillary subsets in trajectory space (Fig. 1e). In contrast, Depp1 is preferentially expressed in pre-Art, consistent with a prior study showing its transient expression in developing arterial endothelial cells and subsequent downregulation in mature vessels 11 (Fig. 1f). Pre-Art and Art express Klf2 and Klf4, genes induced by laminar shear flow and preferentially associated with linear segments compared to branched vessels 9,12 .
Venous EC, comprising HEC and non-HEC (Vn) subsets, share expression of the vein-specifying transcription factor Nr2f2 (Coup-TFII) 13 , and the vein-associated chemokine interceptor Ackr1 (DARC 14 ; Fig. 1f). HEC express genes required for lymphocyte recruitment including Chst4 and Glycam1 15 , with more pronounced expression on HEC distal in tSpace (late HEC; Fig. 1e, f). Consistent with their large size and plump morphology and with prior whole-genome expression studies of sorted HEC 6 , their gene signature is enriched for glycoprotein synthesis (Supplementary Data 1), and they have uniquely high numbers of transcripts per cell (Fig. 1g).
Cells of the Vn subset branch prominently from proximal HEC and TrEC in tSpace projection (Fig. 1e). To identify these cells in the LN vasculature, we imaged whole LN removed shortly after i.v. delivery of fluorescently labeled antibodies (Methods): anti-Ly6c, specific for arteries and capillaries; MECA79 to the Peripheral Node vascular Addressin (PNAd) defining HEV; and anti-PLVAP, which stains capillary and venous EC but not arteries (Fig. 2a). Subset markers allow visualization of arterial entry from the LN hilus, linking to capillary arbors in the cortex which in turn connect to HEC. We identified PNAd − veins downstream of and as a continuation of HEV in the lymph node medulla ( Fig. 2a): they bound injected antibodies to VE-cadherin, PLVAP, and ICAM1 ( Fig. 2b) but were negative for capillary markers Ly6c and podocalyxin (PODXL; Supplementary Fig. 4), as predicted by gene expression (Fig. 2c).
The Vn gene signature includes genes associated with regulation of neutrophil activation (GO:1902563; Fig. 3a) and platelet degranulation (GO:0002576; Supplementary Fig. 3). Vn express Von Willebrand Factor (Vwf), which is stored in Weibel-Palade bodies and released during inflammation to promote platelet adhesion and hemostasis 16 . Surprisingly, Vn lack HEC-associated genes for naïve lymphocyte recruitment (Chst4, Fut7, Ccl21), instead expressing genes for vascular E-and Pselectins (Sele and Selp) and adhesion receptors ICAM1and VCAM1 (Fig. 2c), which mediate myeloid cell recruitment. Neutrophils and monocytes are normally excluded from LN homing, but they enter LN's in large numbers in response to acute inflammation and play an important role in preventing pathogen spread. Prior studies have characterized inflammatory changes in HEC, which enable recruitment of myeloid cells along with lymphocytes [17][18][19] ; but the role of medullary Vn has not been examined. We induced inflammation by footpad injection of the bacterial pathogen S. aureus in mice with green fluorescent protein (GFP) expressing neutrophils (LysM GFP mice). One hour later, mice were transfused with CMTPX-labeled lymphocytes and recruitment was quantified by live two-photon imaging of the popliteal LN. PNAd − medullary venules exhibited massive and exclusive recruitment of GFP + myeloid cells, contrasting with both lymphocyte and induced myeloid cell interactions in HEV (Fig. 3b-e). In contrast to recruitment through HEV which is PNAd-dependent 18 , myeloid recruitment was robust even in the presence of blocking concentrations of anti-PNAd. Anti-Pselectin significantly inhibited GFP + cell accumulation on Vn, and combined inhibition of vascular E-and P-selectins largely abrogated medullary vein interactions of myeloid cells (Fig. 3f), while HEV recruitment was unaffected 18 . Thus venules in the medullary environment have a unique phenotype and function, selectively recruiting myeloid cells to the medulla in acute inflammation. The results illustrate a surprisingly localized BEC programming for differential leukocyte recruitment.  (Fig. 1f). A transitional phenotype subset (TrEC) expresses canonical capillary genes as well as some HEC genes, albeit at low levels compared to bona fide HEC (Fig. 4a, b). TrEC bridge other CapEC to the venous branch in tSpace projections (Fig. 1e), further suggesting a close relation to HEC. They express Chst2 and the HEC genes St3Gal6 and Fut7, which encode glycosyltransferases for the synthesis of sialyl LewisX (SLex) and 6-sulfo-SLex, carbohydrates that can initiate tethering of lymphocytes under shear flow 15 . Synthesis of PNAd, the mature multivalent L-selectin ligand for lymphocyte homing that defines HEV also involves Chst4 and requires the core 2-branching enzyme encoded by Gcnt1. Chst4 and Gcnt1 are nearly undetectable in TrECs, suggesting that TrEC and HEC might display different glycotopes. Thus, we used antibodies to SLex and PNAd to identify TrEC in situ. Imaging revealed a significant population of BEC that co-stained for SLex and for capillary antigens (Fig. 4c, Supplementary Fig. 5) but lacked mature PNAd. They were morphologically thin-walled and were found immediately upstream of HEV, correlating with their position in trajectory space (Fig. 1e). CapEC1 and CapEC2, the two most abundant capillary populations (Fig. 1c), are centrally located among EC in trajectory space, acting as a hub from which differentiated venous and arterial branches emerge (Fig. 1e). The CapEC1 gene signature is enriched for Type IV collagen trimers (GO:0005587; Col4a1, Col4a2) as well as Igfbp3 involved in maintaining vascular tone and blood pressure 20 . They have high vascular endothelial growth factor receptor activity (GO:0005021), yet also express negative regulators of angiogenesis (Hlx and Igfbp3) and have enhanced expression of genes Id1 and Id3 encoding inhibitor of differentiation proteins, bHLH transcriptional regulators that restrain cell differentiation and delay cellular senescence 21,22 (Fig. 1f).
CapEC2 express Egr1, Cxcl1 and genes reflecting NFkB activation (Nfkbia, Nfkbid, Nfkbiz) and JNK activation (Jun, Junb, Jund, Fos, Fosb, Atf3). These genes and pathways are induced by spatial gradients in fluid shear stress [23][24][25] . Compared to CapEC1, CapEC2 also show reduced Gja4 and Id1 and  corresponding to immuno-stained marker proteins; and Sele, Selp and Vcam1 illustrating selective expression by non-HEV vein. Note the decline in Podxl expression from artery to pre-Art to capillary EC subsets, and a corresponding decline in intensity of staining for Ly6c1 as arteries bifurcate into capillaries in situ in a. Mean expression values for each of the four independent cohorts (grey dots) and mean and SEM of the cohort means are also plotted (black diamonds) within the violin plots. enrichment for genes in the MAPK cascade (GO:0000165) and WNT signaling (GO:0060070; Supplementary Fig. 3), characteristics that parallel findings from a ChIP-seq-based analysis of a disturbed oscillatory shear stress response 26 . These CapEC2 features suggest exposure to gradient or disturbed shear stress. A third relatively rare population, CapIfn, has a prominent signature of interferon signaling (GO:0060337), with high expression of transcription factor Irf7, interferon response genes Ifit1, Ifit2, and chemokines CXCL9 and 10 ( Fig. 1f). CapIfn also express Isg15 and Gadd45a, genes associated with an "apoptotic-like" signature described in other single-cell studies 27 .
An angiogenic capillary subset enriched for stem-cellassociated genes. Unsupervised clustering identified a distinctive population of activated capillary endothelial cells, CRP, that occupy their own branch in trajectory space, linking to CapEC1, 2 and TrEC (Fig. 1e). Their gene expression signatures reveal a progenitor-like phenotype. They have uniquely high expression of Ets2 and Sox7, encoding transcription factors implicated in early specification of endothelial cells during development and they display genes and features of stem or progenitor cells in other systems 28 . These developmental EC genes are additionally enriched in early vs late CRP as defined by cell positions along their tSpace branch. CRP also express genes associated with neural and hematopoietic stem or progenitor cells including Cxcr4 29 , Nes 30 , Kit 31 , Lxn 32 , and Sox4 33 ; and they are enriched in spliceosome genes (GO:0097525) including Snrpa1 and Snrpd1 (Fig. 5a), which participate in the acquisition and maintenance of pluripotency in embryonic stem cells 34 . Like embryonic and neural stem cells, CRP uniquely lack expression of Neat1, a long noncoding gene widely expressed in differentiated cells for paraspeckle assembly and double stranded RNA processing during stem cell fate selection 35  signaling than their differentiated progeny, reflecting their developmental plasticity. These characteristics can be quantified by calculating the 'signaling entropy rate' using the SCENT algorithm 36 . Early CRP display higher entropy than other EC subsets, while Art and late HEC subsets display the lowest entropy (Fig. 5b). Dividing cells among CapEC1 share some of these correlates of potency, including high entropy; however, cells with the highest entropy tend to map to the 'origin' of the CRP branch in trajectory space (Fig. 5c). Finally, CRP express Angpt2, Apln, Esm1, Nid2, Pdgfb, and Pgf (Fig. 1f)  that CRPs comprise an activated or primed capillary population with the potential to contribute to vascular maintenance and vessel growth. Consistent with this, CRP are enriched in cells with genes for cell division (Fig. 5d). Although only~10% of CRP have dividing cell signatures, in resting LN~60% of all cells with gene signatures for cell division align with CRP, with most others distributed to CapEC1 and TrEC (Fig. 5e). Interestingly, division increased as Apln expression declined in the transition from early to late CRP and to TrEC, which lack Apln expression (see below).
As CRP selectively express Apln (Fig. 5a) we localized CRP within the LN vasculature of Apln CreER mice (Fig. 5f) by staining for the human estrogen receptor (ER), which serves as a surrogate of Apln expression in these mice 39 . ER + EC were readily visualized as thin-walled endothelial cells in capillaries. HEC were not stained by anti-ER, and importantly ER expression detectable by immunofluorescence histology remained restricted to capillary EC after immune challenge (Fig. 5g, Supplementary Fig. 6). ER + CapEC were rapidly labeled by intravenously injected antibodies to surface antigens, consistent with luminal contact and integration into the capillary endothelium (Fig. 5h). CD276 40 and Cxcr4 38 , known markers of angiogenic potential in EC, are also selectively expressed by CRP (Fig. 5a). Anti-CD276 and CXCR4 antibodies injected intravenously labeled a subset of EC limited to capillaries. In Apln CreER x Rosa26 mTmG (Apln CreER; mTmG ) mice pulsed three days previously with 4-Hydroxytamoxifen (4-OHT) to induce GFP reporter expression in AplnER + cells, GFP + CapEC showed heterogenous staining for injected anti-CD276 and CXCR4 (Fig. 5i).
To assess their fate, we immunized Apln CreER × R26 mTmG mice with Complete Freund's Adjuvant (CFA) 24 h after injection of the short acting tamoxifen metabolite 4-hydroxytamoxifen (4-OHT; serum half-life 6 h 41 ). Three and a half weeks later, many HEC and capillary EC were positive for the reporter, confirming EC neogenesis from AplnERTCre-expressing precursors (Fig. 6). Similar results were obtained in a repeat experiment in which 4-OHT was administered 3 days prior to immunization ( Supplementary Fig. 7) of mice in one leg: GFP expression remained concentrated in scattered CapEC in the control LN even 3 weeks later, whereas numerous GFP + progeny incorporated into HEV in the immunized node. Finally, consistent with the restricted detection of Apln-promoter-driven ERTCre in CapEC even after immunization (Fig. 5f, g), we found that reporter was selectively induced in capillary EC even when 4-OHT was administered 24 h after immunization: Incorporation of reporter + progeny was seen in HEV examined 12 days later ( Supplementary Fig. 8). In conjunction with their expression of genes involved in embryonic vascular development and their enrichment in cells undergoing basal cell division, the results suggest that CRP are a poised regenerative subset or progenitor that contributes to vessel maintenance at steady state and to neogenesis of EC including HEC during LN angiogenesis.
We recently identified an Apln-expressing metaphyseal EC subset in the adult bone marrow that contributes to the hematopoietic stem cell niche 42 : fate mapping in Apln CreER × Rosa26 mTmG mice suggested that the Apln + EC contribute to endothelial homeostasis and neogenesis in the marrow vasculature, as shown here for LN CRP. We therefore sought to identify CRP-like EC in BM and other tissues, taking advantage of scRNAseq profiles of BEC from the Tabula Muris (TM) consortium 43 and the gene expression omnibus (GEO) 44,45 . We performed a global alignment of~38,000 BEC, including cells of our PLN2 and three samples, to cells of our reference PLN1 sample ( Supplementary Fig. 9a). We defined candidate CRP-like cells as cells whose gene expression profiles (1) correlated more highly with the mean gene profile of CRP than with other CapEC or differentiated subsets; and (2) aligned with PLN CRP in UMAP projection. A small percentage (0.3-13%) of BEC met these criteria in most tissues ( Supplementary Fig. 9b): the majority of these expressed Apln as well as other CRP genes including Cd276, Cxcr4, Lxn, Mcam, and Trp53i11 (Supplementary Fig. 9c; Supplementary Data 2). The CRP-like EC express many angiogenic tip-cell-related genes, and like LN CRP share high entropy compared to differentiated BEC in their respective tissues ( Supplementary Fig 9c, d). In addition to aligning with LN capillary CRP, most CRP-like cells express capillary-associated gene markers including Gphibp1, Ly6c1, and Cdh13, and lack venous and artery genes. However, BM CRP-like cells had lower Ly6c1 than cells in other sites. CRP-like EC were extremely rare or absent in liver. In lung, Apln and Ptp4a3 are highly expressed by capillary aerocytes (also known as alveolar type I capillaries 46 and Car4-high EC 47 ). These cells line alveolar spaces, are unique to the pulmonary vasculature and are otherwise unrelated to CRP. The results suggest that CRP are rare but widely distributed EC of capillary phenotype that likely participate in steady state BEC turnover and function as poised progenitors for induced angiogenesis.
Gene regulation along cellular trajectories. We next examined changes in genes and cell features along cellular trajectories. We isolated cells along paths from early CRP to arteries or to HEC, or from CapEC along the venous branch leading to Vn (Fig. 7a), and visualized expression of genes or gene modules by cells along the trajectories (Fig. 7b-e, Supplementary Fig. 10). The heatmap illustrates enrichment of cells with high cell-cycle scores between or along the path from early CRP to CapEC and TrEC, peaking in correspondence with late CRP or the transition to TrEC (Fig. 7b, Supplementary Fig. 10a). Conversely Apln declines rapidly in the transition from early to late CRP and TrEC, and consistent with this dividing CRP have reduced Apln expression compared to early CRP (Fig. 7b, Supplementary Fig. 10a). CapEC on the trajectory to Art express genes implicated in developmental arteriogenesis 48 . These include Sox17, Nrp1, Gata2, Klf2, and genes for Bone Morphogenic Proteins (BMP), Notch and Ephrin signaling components and downstream targets (Acvrl1, Tmem100, Msx1, Notch4 and Notch1, Jag1, Jag2, Hey1, Efnb2) that program the mature arterial phenotype (Fig. 7e). Features reflecting laminar shear stress (Klf2, Klf4, GO:0034616; Fig. 7e) increase along the trajectory to Art, while the oscillatory or gradient shear stress signature peaks in capillary EC, especially the CapEC2-rich region of the trajectory leading to arterial EC (Cxcl1, Egr1, OSS; Fig. 7b, d, e). Genes involved in gas and metabolite exchange (Aqp1 49 , Car2, and Car4) are high in CapEC and lost in the progression to Art. Surprisingly, arterial EC express key genes for elastin fiber assembly (Fbln5, Eln, Ltbp4, and Lox), suggesting that arterial EC may directly participate in the assembly of the inner elastic lamina between EC and smooth muscle cells (Fig. 7d). The CRP to Art trajectory also displays genes involved in arterial EC progenitor migration in development, discussed below.
In contrast, cells along the trajectory from CRP to Vn or HEC express determinants of venous fate and phenotypes. Nr2f2 encodes CoupTF2, a TF required for venous differentiation in development: Nr2f2 is first expressed in TrEC and is maintained along the venous trajectory. Notch is inhibited during venous differentiation 13 , and Notch signaling components are downregulated early along the venous trajectory (Fig. 7d). The lymphotoxin receptor LTBR is required for HEV development and maintenance. Lbtr is broadly expressed, but lymphotoxin beta (Ltb), a component of HEV-specifying lymphotoxin dimers, appears in the late HEC portion of the venous axis, suggesting the potential for autocrine signaling to maintain the HEC phenotype (Fig. 7b). In addition to mapping known mediators and mechanisms of arterial and venular specification to cellular trajectories, the analysis identifies candidate transcription factor genes that may contribute to specialization of large vessels (e.g., Ebf1, Klf9), arteries (e.g., Tsc22d1), HEC (e.g., Xbp1, Meox1), HEC and vein (Aebp1) or veins (Lhx6, Csrnp1, Gata6) (Fig. 7d). For example, Gata6 enhances Tumor Necrosis Factor-alpha induced VCAM1 expression in EC cultures 50 and its expression by TrEC and Vn may thus contribute to selective Vn expression of Vcam1.
G-protein-linked receptors (GPCR) serve as environmental sensors. With their ligands GPCRs regulate vascular development, endothelial function, and migration. We examined expression of GPCR along the aligned trajectories. The vasodilatory flow sensor gene Gpr68 51 first arises in pre-Art capillary EC and is retained in Art (Fig. 7b). Cysltr1, encoding a sensor for myeloid cell-derived leukotrienes 52 , is upregulated in the transition to Vn (Fig. 7b). Genes encoding thrombin and erythropoietin receptors (F2r and Epor), sensors that drive angiogenesis and vasculogenesis 53,54 , are expressed selectively by CRP, suggesting sensitivity to local thrombotic or inflammatory protease activity and oxygen deficiency (Fig. 7b). Patterns of expression of genes for Cxcl12 (in CapEC and pre-Art), and for its chemoattract receptor Cxcr4 (expressed by CRP and subsets of CapEC1 and 2 along the trajectory to Art) suggest retention of a developmental programs for tip-cell migration along capillaries into developing arteries (Fig. 7c), as observed during retinal arteriogenesis 55 . Ackr3 expression by trailing cells is required for Cxcr4/Cxcl12-dependent migration of germ cells in zebrafish development 56 . Ackr3 is expressed by CRP and some CapEC in the arterial branch, suggesting a parallel role. These observations suggest that CRP may be primed for migration. Consistent with this, CRP are enriched in genes for cell locomotion, extracellular matrix remodeling, actin assembly and disassembly and lamellipodium formation (Fig. 7e); but they also express genes for known inhibitors of migration and sprouting behavior including Arhgap18 57 and Csnk2b 58 (Supplementary Data 1).

Discussion
The vascular endothelium plays a central role in lymphoid tissue development and function. Our survey of the LN blood vascular endothelium defines the diversity of EC phenotypes and identifies previously undescribed cell subsets and functions. We uncover multiple subsets of capillary endothelial cells with distinct gene expression including TrEC, a transitional phenotype subset intermediate in gene expression and in physical location between  Fig. 7 Trajectories align EC subsets and mechanisms of EC development and specification to the vasculature. a Cells along KNN-based trajectories were isolated (Methods). b-e Expression of selected genes and gene set enrichment scores along cell trajectories from mature CRP to Art (plotted leftward), and from CRP to HEC or Vn (rightward). Cells along the trajectories were manually gated in the first five principal components of trajectory space and aligned according to distance from early CRP. Representation of cell types along the trajectories is indicated at top. Normalized count data plotted as a function of trajectory distance was smoothed using a gaussian kernel. Trajectory distances were scaled to the longest trajectory (CRP to HEC). Genes were grouped according to biological class or function. Imputed gene expression values for all cells (without duplications) were calculated independently and used for hierarchical clustering within each gene group. Cxcr4, Cxcl12, and Ackr3 are shown as a separate group at the top in c (see results). Average normalized expression values for the min and max subset are indicated to the right of each heatmap. OSS: pooled expression of oscillatory shear stress genes: Egr1, Nfkbia, Junb, Cxcl1. Cell Division: pooled expression of cell-cycle genes 74 . associated with multipotent progenitor cells and participate in basal endothelial proliferation and in vascular neogenesis in response to immunization.
New endothelium in physiologic angiogenesis is thought to arise principally from local EC 59 , and circulating endothelial progenitors do not contribute to new endothelium in immunized LN 4 . Prior studies have highlighted multipotent vessel-resident EC progenitor populations in large vessels 60,61 and lung 62 and precursors in developing bone 63 that can differentiate into diverse EC phenotypes. CRP appear distinct from these populations: They lack Bst1 and Procr, markers of resident endothelial progenitors described in large vessels of liver 60 and fat 61 . They express Kit, a marker of clonally proliferative EC progenitors in lung 64 ; but published single-cell datasets of lung blood vascular EC include only extremely rare EC with CRP-like gene profiles. CRP are also distinct from angiogenic tip cells in that CRP are in direct contact with the lumen and are integral to the vessel lining, whereas tip cells lack luminal contact, instead leading the blind end of invasive sprouts 65 . However, CRP share gene expression and precursor potential with tip cells, and thus may be primed for tip-cell behavior or alternatively for intussusceptive (splitting) angiogenesis. CRP appear closely related to Apln-expressing EC we have observed in the adult bone marrow, which fate mapping suggests contribute to irradiation-induced neogenesis of EC including arterial EC in the marrow compartment 42 . Moreover, we identify rare CRP-like capillary phenotype EC in many tissues. Thus CRP may contribute widely to vascular homeostasis, providing a distributed pool of regenerative cells for local vascular maintenance and replenishment. Consistent with this thesis, a recent scRNAseq study also noted rare Apln + "angiogenic" EC in normal tissues 66 .
Interestingly, in PLN and in most tissues profiled Apln expression appears quite selective for CRP or CRP-like EC. At the gene level, its expression is highest in 'early' CRP, which also show the highest entropy among the BEC profiled here. Apln is downregulated progressively in cells along the trajectory from early to late CRP and CapEC; and its translation to protein as assessed by Apln-promoter-driven expression of Cre-ERT2 remains restricted to capillaries even in immunized lymph nodes, in which many HEV and other BEC are undergoing active proliferation. Thus Apelin is not a general marker of EC activation or proliferation. These considerations emphasize the similarities between early CRP and tip cells in models of sprouting angiogenesis: tip cells like early CRP are Apln high, predominantly non-dividing and display signatures of cell migration; whereas cell division occurs primarily among trailing 'stalk' cells 67 .
The migratory signature of CRP suggests that they may have the capacity to crawl within capillaries to contribute to new vessel formation. We show that CRP express genes for locomotion including a Cxcr4 chemoattractant program with the potential to drive precursor migration toward pre-artery-expressed Cxcl12 for arteriogenesis, as observed for tip-cell progeny in retinal development 68 . We also show CRP expression of Ackr3 encoding the Cxcl12 interceptor Ackr3 (CXCR7): CXCR7 internalizes CXCL12, and in developing cell systems can enhance CXCR4driven directional migration by reducing chemoattractant levels at the "rear" of a migrating population 69 .
CRP display sensory receptors for angiogenesis signals, including thrombin, erythropoietin and VEGFs, and thus appear well programmed to respond to requirements for endothelial proliferation. Consistent with this, we show that nearly half of all BEC with basal cell division signatures align with CRP, predominantly late CRP as just mentioned. CRP may represent a temporary state of activated capillary EC, as suggested by enrichment in dividing cells and similarities to tip cells in their gene expression. Alternatively, they may represent a resident progenitor pool, analogous to regenerative stem or progenitor cells in other settings (i.e., intestinal epithelium and hematopoietic systems). They embody high entropy and express genes characteristic of multipotent progenitors and stem cells. Serial cell transfer or clonal culture will be required to evaluate their clonal progenitor potential and capacity for self-renewal. Our data do not address the precursor potential, in terms of contribution to new vessel formation, of the other EC subsets defined here. CapEC1, TrEC, and early HEC display higher signaling entropy than the most differentiated arterial and venous EC in our samples, suggesting that they may also retain developmental plasticity. Indeed, HEC have been reported to contribute to neosynthesis of HEC and capillary cells when injected into a recipient mouse LN 4 . In the intestinal epithelium recently differentiated cells can even re-acquire multipotency and renew stem cell pools in response to injury or stem cell depletion 70 : Some subsets of BEC may be capable of doing so as well. The relative contributions of different BEC subsets to neosynthesis of specialized EC may be a function of the tissues they reside in and the challenges they experience.
We also observed differences among the major CapEC pool, defining two related clusters, CapEC1 and CapEC2. Both share canonical CapEC genes and genes for gas and metabolite transport, which are downregulated in pre-artery and terminal arterial subsets. CapEC1 feature high expression of inhibitor of DNA binding proteins (Id1, Id3), proteins that act as inhibitors of cell differentiation. CapEC2 express genes (e.g., Egr1, Cxcl1) and pathways (NFkB, JNK, WNT, MAPK cascade) induced by oscillatory or gradient shear stress, which may occur at vessel bifurcations. Interestingly, these genes decline along a trajectory from CapEC2 to arterial EC, while laminar flow associated genes Klf2 and Klf4 progressively increase. Pre-Art may thus correspond to the arterioles and potentially to cells lining arteriovenous communications in the LN.
The ability of medullary veins to selectively recruit myeloid cells, and the selective expression and role of vascular selectins in the process, reveal strikingly local vascular specialization. Previous studies have shown that inflammatory stimuli induce de novo monocyte and neutrophil recruitment to LN but have focused on the role of HEV and the HEV L-selectin ligand PNAd which recruit cells preferentially into the deep cortex (T cell zones). We showed recently for example that neutrophils home via HEV into S. aureus challenged LN in a PNAd-dependent, vascular selectin-independent process 18 . In contrast, we find here that medullary veins lack the machinery for naïve lymphocyte homing, instead selectively recruiting myeloid cells using the vascular selectins. The medulla of the lymph node contains pathogen-trapping lymphatic EC networks that we have recently characterized at the single-cell level 71 : recruitment of neutrophils to the medulla in acute bacterial challenge may contribute to the essential role of myeloid cells in reducing pathogen transit into efferent lymphatics and systemic spread of infection.
Our results show that the major EC subsets, defined by gene signatures, map to specific locations within the steady state vasculature. Indeed, alignment of cells along nearest neighbor trajectories recapitulates the overall architectural arrangement of EC in the blood vasculature. Imaging confirms the computationally predicted positioning of TrEC between CapEC and HEV; of CRP within capillary segments; and the branching of PNAd negative veins from HEV. Trajectory analysis also reveals that capillary EC aligned along trajectories to mature arterial EC express transcriptional programs that, in development, support arteriogenesis. Similarly, genes that program developmental specification of veins are retained along the venous branch. This 'retention' of artery and venous specifying programs may reflect the continuous steady state replenishment and developmental programming of these specialized subsets from dividing CRP or other precursors. Alternatively, retention of developmental genes could serve to pre-program segmental EC differentiation for rapid vascular expansion during immune angiogenesis.
As shown here, single-cell analysis has the potential to identify EC subsets; elucidate developmental processes, transcriptional and regulatory pathways that program their specialization; and map transcriptional phenotypes to the vasculature, providing a molecular blueprint of the vascular endothelium. Targeting specific subsets and processes defined here holds promise to treat a variety of vascular, immune and inflammatory disorders through manipulation of angiogenesis and immune responses. Preparation of lymphoid tissue BECs for flow cytometry. Axillary, inguinal and brachial PLN from 20-30 adult mice: BALB/cJ (PLN1 and PLN2) or mice of mixed background (PLN3) were dissociated as described 6 . To minimize technical variation, in one study (PLN1) male and female PLN were combined before processing the tissue and separated postsequencing using the AddModuleScore function from the Seurat package (v3.1.1) to calculated enrichment of male-specific genes (y-chromosomal genes) and the female-specific gene Xist. Endothelial cells were isolated essentially as described (online version) 6 . In brief, PLN were pooled in HBSS buffer and minced with scissors. Minced tissue was washed 2-3 times with HBSS, resuspended in HBSS media containing 0.2 mg/ml collagenase P, 0.8 mg/ml Dispase II, 0.01 mg/ml DNase, and incubated at 37°C for 10 min. Tissue was then gently disrupted by pipetting up and down several times through P1000 pipette tips with successively smaller bore sizes. Tissue was allowed to settle to the bottom of the tube and cells released by the digestion were transferred to ice cold FBS (final concentration 30%). Fresh digestion buffer was added to the remaining tissue fragments for an additional round of digestion. Digested cells were filtered and hematolymphoid cells were depleted from the cell suspension using anti-CD45 MicroBeads (Miltenyi). Approximately, 5-10 × 10 4 BECs (lin -Gp38 -CD31 + ) were sorted into 100% fetal bovine serum using a FACS Aria (100 µm nozzle;~2500 cells/second). Freshly sorted cell suspensions were diluted with PBS to a final FBS concentration of~10% and centrifuged at 400 × g for 5 min. Supernatant was carefully removed using micropipettes and cell pellets resuspended in the residual volume (~30-50 µl). Cells were counted using a hemocytometer and cell concentration adjusted to 500-1000 cells µl) by addition of PBS with 10% FBS if necessary.
Highly variable genes were identified using the FindVariableGenes function (Seurat, v2.1) 27 . For analyses designed to identify clusters, non-variable genes, cellcycle genes 74 , genes detected in fewer than three cells and genes with an average expression level below 0.3 (normalized and log-transformed counts) were excluded. Supervised cell selection was used to remove cells with non-blood endothelial cell gene signatures: lymphatic endothelial cells (Prox1, Lyve1, Pdpn); Pericytes (Itga7, Pdgfrb); fibroblastic reticular cells (Pdpn, Ccl19, Pdgfra); lymphocytes (Ptprc, Cd52). Top principal components and the FindClusters function (Seurat, v2.1; res = 0.3) were used on a core set of cells (2394) from the PLN1 sample to identify the 8 major clusters. The Arterial, HEC and CRP clusters were further subdivided into Art and Pre-Art, HEC and HEC (late); and CRP and CRP (early), based on canonical marker expression and their position in tSpace projections of PLN1, yielding a total of 11 subsets. The remaining PLN1 cells and cells from the independently processed samples (PLN2 and PLN3) were assigned the identity of the maximally correlated (Pearson) average expression profile of the core PLN1 cell subsets using~3000 common variable genes. Batch effects from technical replicates were removed using the MNN algorithm 75 as implemented in the batchelor package's (v1.0.1) fastMNN function. Cells were classified as dividing or resting using a pooled expression value for cell-cycle genes (Satija Lab Website: regev_lab_cell_cycle_genes). For UMAP and tSpace embeddings, cell-cycle effects were removed by splitting the data into dividing and resting cells and using the fastMNN function to align the dividing cells with their resting counterparts. Dimensionality reduction was performed using the UMAP algorithm and nearest neighbor alignments for trajectory inference and vascular modeling were calculated using the tSpace algorithm 10 . Cells along isolated trajectories were selected by gating within tPC projections 1-5; and illustrated in the figures here. Differential gene expression analysis was performed by comparing each subset to the remaining cells and fitting a zero-inflated negative binomial model using the LineagePulse package, v0.99.20. In order to assess differentiation potency of single cells, the signaling entropy rate (SR) of each cell was calculated 36 after random downsampling of reads to 1000 reads/cell. Data visualization. Heatmaps were generated using the ComplexHeatmap package 76 , scaled to a maximum value of 1. Data for trajectory heatmaps was preprocessed using code adapted from the plot_as_function from cyt (Pe'er Lab): normalized count data was smoothed with respect to trajectory distance using a gaussian kernel and plotted using a discrete color scale. Violin plots were generated using ggplot2; y-axis units for gene expression data correspond to logtransformed normalized counts after imputation. 3d plots were generated using the rgl package v0.100.30, with minor source code modifications for interactive renderings.
GO term analysis. Pooled expression values for GO term gene sets and other sets of genes were calculated as previously described 77 using the AddModule-Score function of Seurat (v3.1.1), which centers values by subtracting pooled expression for random sets of control genes with similar expression levels. To identify biologically relevant GO terms, we first generated a pooled expression matrix by systematically applying the AddModuleScore function to all currently annotated GO terms with at least 3 expressed genes. Differentially regulated GO terms were identified from the resulting "GO term expression matrix" (14300 GO terms by 8832 cells) by comparing each subset to the remaining cells using a Student's t-test. A high degree of overlap was observed with conventional GO term analysis approaches such as analysis of top differentially expressed genes using Enrichr. Because this GO term analysis approach is done on a cell by cell basis, it was particularly useful for the identification of terms whose enrichment spanned multiple subsets, e.g., GO:0015669_gas transport ( Supplementary Fig. 3). Imaging. PLN were imaged following either retroorbital injection of fluorescent labeled antibodies or by fluorescence staining of LN sections. If injected, antibodies (25-75 µg) were administered 5-30 min prior to sacrifice and PLN removal. To image the overall vasculature, the PLN was gently compressed to~35-50 µm thickness on a glass slide. Alternatively, PLNs were fixed with 4% paraformaldehyde, cryoprotected with sucrose, frozen in OCT (Sakura® Finetek) in 2methylbutane (Sigma) on dry ice and stored at −20°C. 50 µm cryo-sections were stained with antibodies according to standard protocols. The slides were imaged using Apotome 2.0 fluorescence microscope or LSM 700 or LSM 880 laser scanning microscope (Zeiss). Noise was removed using Fiji.
For quantification of ER + cells, anti-human estrogen receptor antibody was used as a surrogate stain for Apln in sections from Apln CreER mice. Capillaries, HEV and ER + vessels were enumerated within 50 µm sections at ×20 objective using a grid reticle to determine the relative frequency of each EC subset (1 length unit = 1/8 of the grid height). Sufficient fields were scanned to comprise >5000 HEC assessed for reactivity with anti-ER antibody. Data were expressed as frequency of ER + (as a % of total counted vessels, i.e., length of ER + / length of capillaries + HEC) per lymph node section. We scanned nine LN from two Apln CreER unchallenged mice and nine LN from two Apln CreER mice five days post cutaneous inflammation (one section per LN). The number of HEC per field was determined based on ERG + nuclei within MECA79 + EC using ×10 or ×20 objectives.
Live imaging. LysM GFP mice were injected in the right footpad with 2.5 × 10 7 CFU S. aureus. An hour later mice received intravenously injected Cell Tracker Red CMTPX (Thermofisher) labeled lymphocytes. An hour after lymphocyte injection mice were injected intravenously with a mixture of anti-PNAd − Dylight594 and albumin-Dylight-680 to visualize HEV and vasculature, respectively. Two-photon videomicroscopy was employed to assess the movement and location of neutrophils and lymphocytes in lymph node blood vessels following infection with S. aureus. The right hindlimb popliteal LN was exposed for imaging in the anesthetized mouse. The LN was imaged from 2 to 4 hpi. Image acquisition was performed an upright two-photon microscope (Leica Biosystems TCS SP8 Upright Microscope). For antibody blocking studies, mice were pretreated with blocking antibody 20 min before infection.
Lineage tracing. For lineage tracing of CRP labeled in resting LN prior to immunization, reporter expression was induced in Apln CreER × R26 mTmG mice by i.p. injection of 80 ug/g 4-hydroxytamoxifen (4-OHT; Sigma-Aldrich) either 24 h prior to sacrifice or immunization, or 76 h prior to sacrifice or immunization in separate experiments (see text and figure legends). Immunized mice were sacrificed at the indicated timepoints, and LN taken for imaging. Additional LN to control for leakiness of the reporter included LN from Apln CreER × R26 mTmG mice immunized in parallel without tamoxifen induction. For labeling and tracing of Apln CreERexpressing cells after initial immunization, reporter expression was induced in Apln CreER × R26 tdTomato mice the day after oxazolone skin painting, and lymph nodes were imaged 24 h later or after 11 days.
Statistical analysis. Statistical significance between two groups was calculated using a two-way ANOVA corrected with Tukey. The rule of three was applied to determine 95% confidence intervals for the enumeration of ER + cells. A likelihood ratio test was used for differential gene expression analysis assuming an underlying zero-inflated negative binomial distribution. P-values were adjusted for multiple comparisons by calculating the false discovery rate (FDR) and adjusted p-values < 0.001 were considered significant.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The authors declare that all data supporting the findings of this study are available within the article and its Supplementary Information files or from the corresponding author upon reasonable request. Raw single-cell RNAseq data have been deposited in the Gene Expression Ominibus (GEO) database under accession code: GSE140348. Processed single-cell data used in Supplementary Fig. 9 were obtained from the GEO (GSE108892 bone marrow endothelial cells; GSE106514 peripheral lymph node high endothelial cells) and the Tabula Muris consortium [https://registry.opendata.aws/tabula-muris-senis/]. Interactive web apps to visualize single-cell RNAseq data are available the PLN1 dataset [https://stanford.io/2qzJ8Hl] and the integrated dataset (PLN1, PLN2, and PLN3) [https://stanford.io/2WXR811].

Code availability
Code used in this study, including a customized implementation of the MAGIC algorithm, are available as part of the magicBatch package [https://github.com/kbrulois/ magicBatch].