A single-cell transcriptional roadmap of the mouse and human lymph node lymphatic vasculature

Single-cell transcriptomics promises to revolutionize our understanding of the vasculature. Emerging computational methods applied to high dimensional single cell data allow integration of results between samples and species, and illuminate the diversity and underlying developmental and architectural organization of cell populations. Here, we illustrate these methods in analysis of mouse lymph node (LN) lymphatic endothelial cells (LEC) at single cell resolution. Clustering identifies five well-delineated subsets, including two medullary sinus subsets not recognized previously as distinct. Nearest neighbor alignments in trajectory space position the major subsets in a sequence that recapitulates known and suggests novel features of LN lymphatic organization, providing a transcriptional map of the lymphatic endothelial niches and of the transitions between them. Differences in gene expression reveal specialized programs for (1) subcapsular ceiling endothelial interactions with the capsule connective tissue and cells, (2) subcapsular floor regulation of lymph borne cell entry into the LN parenchyma and antigen presentation, and (3) medullary subset specialization for pathogen interactions and LN remodeling. LEC of the subcapsular sinus floor and medulla, which represent major sites of cell entry and exit from the LN parenchyma respectively, respond robustly to oxazolone inflammation challenge with enriched signaling pathways that converge on both innate and adaptive immune responses. Integration of mouse and human single-cell profiles reveals a conserved cross-species pattern of lymphatic vascular niches and gene expression, as well as specialized human subsets and genes unique to each species. The examples provided demonstrate the power of single-cell analysis in elucidating endothelial cell heterogeneity, vascular organization and endothelial cell responses. We discuss the findings from the perspective of LEC functions in relation to niche formations in the unique stromal and highly immunological environment of the LN. Highlights Computational alignments (“trajectories”) predict LN LEC organization in situ, revealing a continuum of phenotypes punctuated by specialized clusters Multiple intermediate phenotypes suggest LEC malleability Gene profiles define niche-specific functional specialization Medullary sinus LECs are comprised of Ptx3-LECs and Marco-LECs Distinct mechanisms for pathogen interactions and matrix modeling Ptx3-LECs: paracortical and central medullary sinuses near hilus; enriched for genes driving lymphangiogenic responses and lymphocyte egress Marco-LECs: peri-follicular medullary sinuses; macrophage-associated genes, complement and coagulation cascade Niche-specific responses to inflammation IFN gene responses in SCS floor and medullary sinus LECs Suppression of LEC identity genes in responding subsets Conserved and unique LEC subsets and gene programs across species Core subsets common to mouse and human Greater diversity of subsets and intermediates in human LN LECs

Single-cell transcriptomics promise to revolutionize our understanding of the vasculature. Emerging computational methods applied to high-dimensional single-cell data allow integration of results between samples and species and illuminate the diversity and underlying developmental and architectural organization of cell populations. Here, we illustrate these methods in the analysis of mouse lymph node (LN) lymphatic endothelial cells (LEC) at single-cell resolution. Clustering identifies five well-delineated subsets, including two medullary sinus subsets not previously recognized as distinct. Nearest neighbor alignments in trajectory space position the major subsets in a sequence that recapitulates the known features and suggests novel features of LN lymphatic organization, providing a transcriptional map of the lymphatic endothelial niches and of the transitions between them. Differences in gene expression reveal specialized programs for (1) subcapsular ceiling endothelial interactions with the capsule connective tissue and cells; (2) subcapsular floor regulation of lymph borne cell entry into the LN parenchyma and antigen presentation; and (3) pathogen interactions and (4) LN remodeling in distinct medullary subsets. LEC of the subcapsular sinus floor and medulla, which represent major sites of cell entry and exit from the LN parenchyma respectively, respond robustly to oxazolone inflammation challenge with enriched signaling pathways that converge on both innate and adaptive immune responses. Integration of mouse and human single-cell profiles reveals a conserved cross-species pattern of lymphatic vascular niches and gene expression, as well as specialized human subsets and genes unique to each species. The examples provided demonstrate the power of single-cell analysis in elucidating endothelial cell heterogeneity, vascular organization, and endothelial cell responses. We discuss the findings from the perspective of LEC functions in relation to niche formations in the unique stromal and highly immunological environment of the LN.

INTRODUCTION
Lymph nodes (LNs) serve as hubs for the interaction and communication between tissue-derived and blood-derived immune cells (1). Integrated along the large collecting lymphatic vessels, they are vital sensors of tissue damage, constantly sampling the incoming lymph (2). The LN comprises a complex network of lymphatic sinuses surrounding a dense parenchyma, which mainly consists of immune cells but also specialized blood vessels and a network of mesenchymal cells (1)(2)(3). Segregated B cell (cortex) and T cell (paracortex) areas characterize the LN architecture (4). It is well-established that the LN stromal cells play a central role in maintaining both this structure and the immunological functions of the LN, providing chemotactic cues, cytokines, and a structural reticular framework that guide immune cell positioning, migration, survival, and activation [reviews: (3,5)]. Single-cell sequencing has enabled delineation of nine distinct clusters of murine LN mesenchymal cell phenotypes (6) underlining the complexity needed to maintain the LN structure and coordinate immunity.
The lymphatic vasculature is the first structural component of the LN encountered by incoming lymph-borne molecules or cells. Recent studies have revealed an intriguing regional specialization and cellular heterogeneity that characterize the LN lymphatic endothelium and differentiate the LN lymphatic endothelial cells (LECs) from LECs in peripheral lymphatic vessels (7)(8)(9)(10)(11)(12). Subset-specific markers with functional implications include the atypical chemokine receptor ACKR4 (also known as CCRL1), specifically expressed by the LEC layer that forms the ceiling (cLECs) of the subcapsular sinus (SCS) (9), where lymph enters from the afferent collecting vessels. ACKR4 is a scavenger receptor for the homeostatic chemokines CCL19, CCL21, and CCL25 (13) and controls entry of tissue-derived dendritic cells (DCs) into the LN through controlling the formation of CCL21 chemokine gradients across the SCS (9). Leukocyte entry occurs primarily through the SCS floor LECs (fLECs), which in mouse express the mucosal vascular addressin cell adhesion molecule 1 (MAdCAM-1) (7,14) among other adhesion and attractant molecules that can control leukocyte transmigration (7,15). The SCS also functions as a physical barrier and gateway, enabling size-restricted access of antigens to the LN parenchyma (16): the glycoprotein plasmalemma vesicle-associated protein (PLVAP), together with cLEC-expressed caveolin 1 (CAV1) (10), form sieve-like diaphragms in the transendothelial channels that bridge the SCS to the conduit system and descend from the SCS floor (10). This structural barrier is complemented by a dense network of macrophages closely associated with SCS, providing essential innate immune functions and a filtering system for pathogens (17). The SCS and medullary sinus macrophage niches in the LN were recently shown to be dependent on LEC-expressed CSF-1 (colony stimulating factor 1) (18) and receptor activator of nuclear factor κ B (RANK)/RANKL signaling between the LN LECs and SCS lining mesenchymal cells (19).
LN LECs can also directly influence adaptive immune responses, either through presentation of tissue antigens, which contributes to the maintenance of peripheral tolerance (20)(21)(22), or by serving as reservoirs for antigens (23). LN LECs express the immune-checkpoint ligand programmed death-ligand 1 (PD-L1) (also known as CD274) (7,24), an inhibitor of T cell activation, and lack expression of co-stimulatory genes (24), which may explain their role in tolerance. PD-L1 is expressed selectively in the floor of the SCS (fLECs), cortical sinuses, and parts of the medulla (7,24). Genes that influence the communication between LECs and their surroundings could contribute to endothelial regulation as well, and interestingly, PD-L1 expression moderates proliferation and enhances survival of LN LECs in inflammation (25). The diverse and site-specific specialization of the LN lymphatic endothelium is at least partly dependent on crosstalk with immune cells, with contributions from B cells, T cells (7,26), and mesenchymal stromal cells (19). Hence, the LNs provide both a unique model system to explore endothelial cell interactions with their surroundings and a model for exploring endothelial diversity and phenotypic plasticity.
Our recent single-cell analysis of the human LN LECs revealed the complexity and specialization of the LN lymphatic endothelium in man (27). A detailed profiling of the mouse lymphatic endothelium and comparison of human and mouse LN lymphatic endothelium is still missing. Here we provide single-cell transcriptomic analysis of the mouse LN lymphatic vasculature. We show that computational alignments (relationships revealed by nearest neighbor trajectory inference) recapitulate key aspects of the tissue architecture and predict physical relationships between LN LECs in the tissue, illustrating the power of single-cell analysis for understanding the organization of the vascular endothelium. Cross-species analyses further allowed us to define conserved and divergent LEC phenotypes and lymphatic vascular niches. Notably, the analysis delineates two specialized subsets of medullary sinus LECs, present in both mouse and human that are distinct in gene expression and location in the LN. Using the mouse LN response to cutaneous oxazolone as a model, we show that sites of immune cell entry and exit from the LN, fLECs, and paracortical sinuses, respectively, respond rapidly to inflammation, whereas structurally important cLECs are less affected. Together, the results demonstrate the power of bioinformatic tools for elucidating endothelial cell heterogeneity, physical relationships, and cellular responses in complex vascular beds, and provide a basis for future detailed analysis of human and mouse LN LEC responses in disease.

Single-Cell Trajectories Model Tissue Architecture and Physical Relationships Between LN LECs
To assess the heterogeneity of the lymphatic vasculature and the relationships between lymphatic vascular niches, we analyzed lin − podoplanin (PDPN) + CD31 + cells from mouse peripheral LNs (i.e., axillary, brachial and inguinal) by singlecell RNA sequencing (scRNA-seq) using the 10x Genomics system ( Figure 1A). An independent set of LN LECs was sorted as single cells and subjected to SMART-seq2 analysis (28) (Figure 1A). Prospero homeobox protein 1 (Prox1) was used as a pan-LEC marker in the analysis (29). Blood endothelial cells (Flt1 + Ly6c1 + ), immune cells (Cd45 + ), and pericytes (Pdgfra + or Pdgfrb + ) were identified by marker gene expression and excluded from further analysis. We used a combination of unsupervised clustering and graph-based methods to determine LEC subsets ( Figure 1B; Methods).
In addition to gene profile-based dimensionality reduction with t-distributed Stochastic Neighbor Embedding (tSNE) ( Figure 1B) and Uniform Manifold Approximation and Projection (UMAP) (see below), we used a trajectory analysis algorithm, tSpace (30), to define high dimensional nearest neighbor alignments that emphasize the continuum of cell phenotypes and reveal transitions between related cells ( Figure 1C). In developing cell systems, trajectory inference methods can model developmental sequence ("pseudotime") (30). In the resting adult vascular endothelium, they facilitate computational modeling of vessel architecture (27,31). Subsets and alignments were shared in mice with different genetic backgrounds (C57BL/6 and BALB/c) ( Figure S1A). Cells from both male and female mice were represented in the dataset and comprised similar phenotypes ( Figure S1B). However, the study was not designed or powered to examine strain-or sex-dependent differences in LN LECs, and these remain to be further explored.
Our analyses identified major clusters representing LN cLECs and fLECs based on the expression of known SCS ceiling marker Ackr4 (9) as well as SCS floor markers Madcam1 (7,14) and the chemokine Ccl20 (15,27) (Figure 1D and Figure S1C). Candidate valve-related LECs display high expression of known lymphatic valve markers, including the transcription factor Forkhead box protein C2 (Foxc2) (32) (Figure 1D and Figure S1C). A minor cluster, most apparent in tSpace projections ( Figure 1C) and associated with valve, expressed higher levels of peripheral lymphatic vessels markers including lymphatic vessel endothelial hyaluronan receptor 1 (Lyve1) and the chemokine Ccl21, together with lower expression of Foxc2 compared to the candidate valve LECs (Figure 1D). This population likely represents collector or pre-collector LECs (2,32). These subsets were poorly represented in our data and we discuss them in the context of the architecture of LN LECs, but exclude them from detailed differential gene expression analyses below. Unsupervised clusters separate Lyve1 high candidate medullary sinus LEC into two subsets, referred to here as Ptx3-LECs (Ptx3 + ) and Marco-LECs (Marco + ) ( Figure 1D and Figure S1C).
Alignment of LECs in tSpace visualization recapitulates known connections within the complex lymphatic endothelial network in LN and reveals previously unappreciated relationships as well as intermediate, putative transitional phenotypes. In addition to the subsets identified by clustering, we highlight two transitional populations here, for discussion below. Transition zone tzLECs comprise a link between fLECs and Marco-LECs in tSpace ( Figure 1C). Additionally, a minor "bridge" population (B) aligns along a direct path between fLECs and cLECs in tSpace projection ( Figure 1C). As emphasized by numbering of subsets in Figure 1C and illustrated schematically in Figure 5A, trajectories starting from candidate collecting LECs (1) lead prominently to valve (2) and to SCS ceiling LEC (3), consistent with their known physical connections. cLECs branch to fLECs through the bridging population (B) but also transition prominently to Ptx3-LECs (4), Ptx3-LECs to Marco-LECs (5), and Marco-LECs via tzLECs (6) to fLEC (7) along a phenotypic sequence or path that is well-represented in all LEC datasets here.

In situ Localization of LN LEC Subsets
To define LEC subsets and their niches in situ, we carried out immunofluorescence staining using antibodies to subset differentiating markers predicted from gene expression ( Figure 1D) in inguinal LNs of Prox1-GFP reporter mice (33), where all LECs express the fluorescent reporter GFP (Figure 2). scRNA-seq predicts that LYVE-1 is absent in cLECs and low in fLECs, but is highly expressed in Marco-LECs and Ptx3-LECs ( Figure 1D). Staining for LYVE-1 highlights the cortical, paracortical, and entire medullary region of mouse inguinal LN (Figures 2A,B). Expression of Cd274 (Figure 1D), which encodes the immune checkpoint inhibitor PD-L1, distinguishes fLECs, Marco-LECs, and tzLECs (Cd274 high ) from Ptx3-LECs (Cd274 low ) ( Figure 1D and Figure S1C). Co-immunostaining of PD-L1 and LYVE-1 shows that PD-L1 expression defines the fLECs and discrete regions of LYVE-1 high medullary sinus LECs that are distant from the hilus (Figure 2A). Expression of Marco, which encodes the scavenger receptor Macrophage Receptor with Collagenous Structure (MARCO), is highly selective for Marco-LECs ( Figure 1D): MARCO is a known medullary sinus marker in human LN (27,34). We show that MARCO expression pattern in the mouse LN medulla selectively overlaps with that of PD-L1 ( Figure 2B). Conversely, in situ hybridization for Ptx3 (RNAscope) specifically outlines the central medullary sinuses with particularly strong signals in the paracortical cords ( Figures 3A,B), overlapping with the area of LYVE-1 high PD-L1 − MARCO − cords outlined by protein staining (Figure 2). Co-staining of MARCO and the B-cell marker B220 shows that MARCO high regions start from the cortical areas (defined based  predicts that Ptx3-LECs are contiguous with the SCS ceiling ( Figure 1C), which is also evident along the efferent side of the LN where the LYVE-1 − outermost endothelial layers (cLECs) connect to the LYVE-1 + cords (i.e., Ptx3-LECs) ( Figure S3).
As noted above, trajectory analysis predicts a close relationship between Marco-LECs and fLECs, and identifies a transitional population tzLECs between them. tzLECs are characterized by variable and often intermediate expression of fLEC (e.g., Madcam1 and Ccl20) and Marco-LEC associated genes (e.g. Marco and Lyve1) ( Figure 1D). Histologic analysis identifies a region of LECs between fLECs and Marco-LECs that displays expression patterns consistent with this transitional phenotype. Cells of this phenotype can be observed between the fLECs and Marco-LECs in the region between the bilateral lobes of the inguinal LN and at the LN margins, where the lining cells display reduced MAdCAM-1 compared to the majority of fLEC and increasing levels of LYVE-1 (Figure 4). Few if any genes are specific to tzLECs (i.e., not shared with fLECs or Marco-LECs), and unbiased non-negative matrix factorization does not identify a gene set specific to the subset (data not shown). Thus, the tzLECs highlighted here are part of a continuum between fLEC and Marco-LEC populations, defining an area of vascular zonation. Consistent with computational predictions, the results favor a model in which Marco-LECs occupy the perifollicular medulla in association to the LN cortex, with the transitional phenotype tzLECs bridging them to the SCS floor ( Figure 5A).
Hence, we have demonstrated the organization of LEC in situ, with SCS ceiling converging on Ptx3-LECs, Ptx3-LECs to Marco-LECs and tzLECs, then leading to fLECs (Figures 2-4), which mirrors the phenotypic progression revealed by trajectory analysis (Figure 1C). A schematic of the LN LEC connections is shown in Figure 5A.

LEC Molecular Phenotypes Correspond With Their Vascular Niches and Functions
The cLEC Molecular Profile Reflects Interaction With the LN Capsule cLECs connect to afferent and efferent collecting vessels, and coexpress several genes with valve-related LECs, including higher levels of Foxc2 (36) compared to other LN subsets ( Figure 1D and Figure S1C). Deletion of Foxc2 in lymphatic vessels has been shown to lead to downregulation of the cLEC marker Ackr4, but also to mislocation of surrounding smooth muscle cells (36). cLECs interact physically with the relatively rigid capsule that surrounds LNs. The capsule is a dense network of connective tissue and hence is a unique extracellular matrix (ECM), distinct from other parts of the LN (4). Reflecting this, cLECs are highly enriched for genes encoding ECM proteins, Multimerin protein 1 and 2 (Mmrn1 and Mmrn2), together with high levels of Platelet-derived growth factor a and b (Pdgfa and Pdgfb), Jagged-1 (Jag1) (a ligand for Notch receptors on both endothelial and smooth muscle cells), and Endothelin 1 (Edn1) (a vasoactive peptide and known lymphatic constrictor), proteins known to maintain interaction between endothelial cells and mural cells (37,38) (Figure 5B and Figure S4). The capsule was recently shown to have a CD34 + stromal cell subset with expression of CD248 (6), a ligand for MMRN2. MMRN2 can also interact with cLEC expressed MMRN2-ligand Cd93, which was recently coupled to ECM fibrillogenesis in tumor blood endothelium (39). cLECs express the BMP family members Growth differentiation factor 10 (Gdf10) and the Bone morphogenic protein 4 (Bmp4), which differentiate cLECs from fLECs as well as bridging cells in the SCS ( Figure 5B and Figures S4, S5). fLEC instead express Bmp2 ( Figure 5B and Figures S4, S5). The expression of BMPs may mediate both autocrine and paracrine signaling to the surrounding stromal and immune cells.   Figure S4). These adhesion receptors may also help retain the closely associated SCS macrophages (17). fLECs are also enriched for chemokines, including the known fLEC marker Ccl20 (CCR6 ligand) (15,27) but also Cxcl1 (CXCR2 ligand) ( Figures 1D, 5B and Figure S4). CCL20 has been linked to innate lymphoid cell (ILC) trafficking across the SCS (15) and may affect the cross-talk between CCR6-positive DCs and LECs in antigen presentation (40). CCL20 could also influence B cell homeostasis, as memory B cell precursors are distinguished by high CCR6 expression in lymphoid organs of both mouse and man (41) and are closely associated to the SCS (42). CXCL1 may instead contribute to neutrophil migration and recruitment (43). fLECs and Marco-LECs, both macrophage-rich niches (17), share high expression of Csf1 ( Figure 5B and Figure S4), and fLECs express Platelet factor 4 (Pf4) (also known as CXCL4) ( Figure 5B and Figure S4). CSF-1 and CXCL4 both regulate monocyte and macrophage functions (18,44) and LEC-derived CSF-1 is required to maintain the LN macrophage niches (18). As mentioned earlier, Bmp2 is specifically expressed by fLECs ( Figure 5B and Figures S4, S5) that are also enriched for both Smad1 and Smad4 ( Figure 5B and Figure S4) suggesting a role for BMP-signaling in fLEC homeostasis.
Another distinguishing feature of fLECs is their expression of CD74 ( Figure 5B and Figure S4), the MHC class IIassociated invariant chain. CD74 is involved in the formation and transport of MHC class II-antigen complexes (45). Together with high levels of H2-Ab1 ( Figure 4B and Figure S4), this supports a major role of this subset in LN LEC-mediated antigen presentation. Together with the high expression of PD-L1 (Cd274) in fLECs (Figures 2, 5B and Figure S4), and given that antigens are continuously transported from afferent lymphatics into the SCS, this places the fLECs as the major niche for LN LEC-mediated tolerance.
In contrast to all other subsets, fLECs and tzLECs lacks detectable expression of the scavenger receptor Stabilin 2 (Stab2), a sinusoidal endothelial marker (46), and have low expression of the tight junction protein Claudin 5 (Cldn5) ( Figure 5B and Figures S4, S6). The latter suggests that junctional properties of the fLECs differ compared to other LN LECs and may reflect requirements for the active immune migration at this site.

Marco-LECs: MARCO + Medullary Sinus LECs
Immunofluorescence staining for MARCO in mouse inguinal LN shows that it delineates the medullary sinuses adjacent to the B cell follicles in the LN cortex (perifollicular sinuses). MARCO expression coincides with PD-L1 expression in the medullary sinuses, and is excluded from the PD-L1 negative medullary sinuses closest to the hilus (Figure 2 and Figure S2). Lymph from SCS passes through this lymphatic meshwork, abundant in macrophages that can capture lymph-borne pathogens (17,47). While MARCO is also marker for medullary macrophages (47), a major part of the MARCO expression within the adult LN is however confined to LECs (Figure 2). MARCO expression is also abundant in human LN medullary sinuses (27,34). In macrophages, MARCO expression facilitates phagocytic clearance by binding both gram-negative and positive bacteria (48), but its function in LECs is not known. Marco-LECs also express the C-type lectin Clec1b (CLEC2) ( Figure 5B and Figure S4), a ligand for podoplanin (PDPN). CLEC2 is also expressed by myeloid cells, including DCs (49), and interaction between PDPN + fibroblastic reticular cells (FRCs) and CLEC2 + DCs has been implicated in FRC contractility and LN expansion in inflammation (49,50). LN LEC-expressed CLEC2 may mediate both homotypic and heterotypic cell interactions, since PDPN is highly expressed by LECs and by surrounding FRCs. Stab2 is highly expressed by Marco-LECs and Ptx3-LECs ( Figure S6), and together with Ptx3-LECs they are also enriched for the coagulation factor VIII (F8) ( Figure 5B and Figure S4). We previously showed that LN LECs are a major source of FVIII production (51). The present result shows that F8 along with Stab2 are features of medullary LEC subsets. Their co-expression may be a general feature of sinusoidal EC, since both are also characteristic of hepatic sinusoidal endothelium. FVIII forms a complex with von Willebrand factor (vWF) in plasma, a complex that is cleared by STAB2-expressing liver sinusoidal blood endothelial cells in a vWF-dependent manner (52). However, the very low expression of vWF by mouse LN LECs (51) (and data not shown) suggests that FVIII may function as a procoagulant factor independent of vWF, possibly promoting the formation of fibrin mesh to block the spread of pathogens such as Staphylococcus aureus in the medulla (53).
The results show that Marco-LECs share gene expression patterns with myeloid cells, suggesting that this LEC subpopulation has a major role in innate immune functions. Interestingly, Marco-LECs and Ptx3-LECs express higher levels of Kdr, encoding VEGFR2, than other LN subsets ( Figure 5B and Figure S4), which is discussed in further detail below.
Ptx3-LECs share gene expression with dermal capillary LECs, with which they share morphologic features. Both have sproutlike blind ends specialized for attracting leukocytes and fluid flow into the lymphatics (4, 55) (Figures 2, 3). Ptx3 itself is not expressed by dermal capillary LEC in mouse, but capillary and Ptx3-LECs share Itih5, Mrc1, Ccl21, and Lyve1 (56-58) (and data not shown). CCL21 mediates recruitment of CCR7 + dendritic cells (DCs) and T cells into dermal lymphatic capillaries (58,59). Although in the peripheral LN Ptx3-LECs analyzed here Ccl21 expression is relatively low (compared with levels in candidate collecting vessels in Figure 1D), it is more highly expressed in Ptx3-LEC in mesenteric LN (not shown) and in the human LN, as discussed below. Thus, CCL21 in Ptx3-LECs may contribute to leukocyte-LEC interactions and egress of CCR7 + cells including naïve B and T cells, possibly complementing the known role of sphingosine-1-phosphate (S1P) in this process (55,60). Expression of the hyaluronan receptor LYVE-1 promotes the transmigration of DCs into blind-ended capillary lymphatic vessels (57). MRC-1, a known marker for subsets of macrophages, binds microorganisms and promotes phagocytosis (48). MRC-1 expression in peripheral lymphatic vessels also facilitates interaction with CD44-expressing immune cells (56), and could thus interact with activated CD44 high T cells in egress from the LN. Thus, the paracortical and medullary sinus Ptx3-LECs share with peripheral capillary LECs gene programs for LEC-immune cell interactions, and support lymphocyte entry into medullary sinuses and exit from the LN.
Ptx3 itself, as well as other genes selective to Ptx3-LECs, have known or proposed roles in the extracellular matrix (ECM) (54). PTX3 binds collagen and fibrinogen-domain containing proteins, including ECM components, but also other pattern recognition molecules like C1q and ficolins (61); thus, it may amplify innate pathogen recognition mechanisms and contribute to maintenance and repair of the LN medullary environment. The Ptx3-LEC marker Reelin (Reln) is an extracellular matrix glycoprotein. In the periphery, Reelin is associated with the transition of peripheral capillary vessels to collecting vessels and contributes to the cross-talk between collecting vessel and surrounding smooth-muscle cells (62). ITIH5 was originally isolated in a complex with hyaluronan (63) which can stabilize the ECM (64). Notably, Ptx3-LECs also express high levels of the hyaluronan receptor LYVE-1 and ITIH-proteins have been shown to interact with PTX3 (65). This provides a possible functional link between LYVE-1, ITIH5, and PTX3 in the Ptx3-LEC niche. Ptx3-LECs are also enriched for the collagen and integrin-binding ECM protein TGFbeta-induced protein (Tgfbi) (Figure 5B and Figure S4), which is induced in LECs by hypoxia (66). Together with high levels of the ECM protein Mmrn1 (Figure 5B and Figure S4), these genes support the notion of a highly specialized medullary ECM.
Ptx3-LECs show the highest levels of Vascular endothelial growth factor receptor 3 (Vegfr3, or Flt4) (67) and of its coreceptor Neuropilin 2 (Nrp2) (68) (Figure 5B and Figure S4), suggesting a higher responsiveness to vascular endothelial growth factor C (VEGFC) (69,70). VEGFC is a major driver of lymphangiogenesis in development (69) and in the adult lymphatic vasculature (71). VEGFR2 (Kdr), which binds both VEGFC and VEGFA (67), is also highly expressed in Ptx3-LECs, although it is most enriched in Marco-LECs. Thus, VEGFA and VEGFC may elicit niche-specific responses in medullary sinus LECs, acting differently on Vegfr3 high Vegfr2 intermed Ptx3-LECs and Vegfr3 intermed/low Vegfr2 high Marco-LECs. Interestingly, while SCS fLEC and cLECs are established in embryogenesis (36), the medullary sinuses form by sprouting during the first two postnatal weeks (72). Thus, high Vegfr3 (Flt4) in Ptx3-LEC could reflect retention of their postnatal programs for VEGFC-dependent sprouting (72), and may imbue them with the capacity to regenerate or expand rapidly in response to the requirements of LN inflammatory responses. It is wellestablished that it is primarily the LYVE-1 + medullary and paracortical sinuses that expand in response to inflammation and in tumor draining LNs (73,74). The Ptx3-LEC marker Angiopoietin 2 (Angpt2) has also been linked to LEC proliferation (75,76): it can induce lymphatic hyperplasia when overexpressed (76) and is highly expressed by newly formed lymphatic vessels in inflammation-induced lymphangiogenesis (75). Ptx3-LECs, like cLECs, but in contrast to all other parenchymal LN LEC subsets, also lack expression of PD-L1 (Cd274). Loss of PD-L1 is associated with increased LN LEC proliferation but also increased LEC apoptosis after inflammation (25). Thus, low PD-L1 suggest higher responsiveness not only to proliferationinducing signals but also higher sensitivity to apoptotic cell death, the latter thought to facilitate LN contraction after inflammatory stimuli (25). Taken together, high Vegfr3, Vegfr2, Nrp2, Angpt2, and low Cd274 (PD-L1) suggest a major role of Ptx3-LECs in LN remodeling.

A Cellular Bridge Between the SCS Ceiling and Floor
As discussed earlier, trajectory analysis not only predicts a transitional population, tzLECs, that link fLEC and Marco-LEC, but further identifies a bridge population that connects cLEC and fLEC. This bridge population is variably represented in the mouse samples ( Figure 1B and Figure S1A) but is prominent in the human (see below, Figure 10). We plotted gene expression of cells along the path from cLEC to fLEC (Figure S5 and Figure S7A) and found that individual cells within the bridge display variable levels of cLEC and fLEC marker transcripts including Lyve1 (Figure S7A), but overall a decreasing gradient of cLEC markers (e.g., Ackr4 and Cav1) and an increasing gradient of fLEC markers (Madcam1, Pd-l1/Cd274) in the progression from cLEC to fLEC ( Figure S7A). We speculate that these cells participate in formation of previously described "transendothelial channels, " physical bridges that link the subcapsular ceiling and floor (10,12,36). Immunofluorescence imaging indicate that these cords consist not only of channels (10,36), but also nucleated LECs traversing the sinus, as indicated by nuclear staining for the lymphatic reporter Prox1-GFP ( Figure S7B). We used immunofluorescence and in situ hybridization for RNA (RNAscope) to evaluate bridge cell expression of cLEC and fLEC markers ( Figure S7 and Figure S5). Most bridging cells displayed defining phenotypes of cLEC (LYVE-1 − , PD-L1 − , MAdCAM-1 − ) or fLEC (LYVE-1 + , PD-L1 + , MAdCAM-1 + ) on the protein level. Cells co-expressing these cLEC and fLEC markers (as predicted for cells of the computationally defined bridge) were rare or absent ( Figure S7). However, RNAscope revealed that bridging cells lack mRNA for Bmp4 (cLEC marker) and have intermediate expression of Bmp2 (fLEC marker), a pattern consistent with the single cell profiles of the predicted bridge ( Figure S5). Although this observation provides some support for the correspondence between the physical and trajectory-inferred cellular bridges, they are far from being conclusive. An alternative possibility is that fLEC and cLEC can interconvert in response to local environmental signals, and that the transcription profiles of these "bridge" cells are snapshots of cells in process of these transitions. Further experiments are required to map the computationally defined "bridge" cells within the LN lymphatic network and to confidently identify the profiles of the in situ cellular links between SCS floor and ceiling.

Niche-Specific Inflammatory Responses in LN LECs
Oxazolone (Oxa) is a chemical compound used in a prototypical model for inflammation-induced LN hypertrophy (73). We here chose to look at initial LEC responses 48 h after topical skin application of Oxa, before any major lymphangiogenic responses are expected. Lymphangiogenesis is a late effect in LN remodeling, peaking at 5-7 days after Oxa, or other inflammatory stimuli (25,73,77). We aligned LEC of the Oxa-treated mouse with the untreated control of the same genetic background and observed two new subsets, termed Ox1 and Ox2, in the immunized group ( Figure 6A).
Inspection of the tSpace projections tUMAP (UMAP on trajectory distances) and tPC (principal components of trajectory distances) suggests that Ox1 arises from fLECs, and Ox2 from Ptx3-LECs ( Figure 6A). Pairwise Pearson correlation based on the top 1,000 variable genes confirmed that Ox1 is most correlated with fLECs (Pearson correlation coefficient = 0.97) and Ox2 with Ptx3-LECs (Pearson correlation coefficient = 0.92) (Figure 6B). We asked which genes were differentially expressed in Ox1 and Ox2 due to Oxa treatment ( Figure 5C). Compared to fLEC and Ptx3-LEC in untreated LN, respectively, we found that both Ox1 and Ox2 displayed typical interferon (IFN) responses, with dramatic induction of Cxcl9 and Cxcl10 as well as Irf7, the master regulator of type-I IFN-dependent immune responses (78) (Figure 6C). CXCL9 and CXCL10 bind CXCR3, which mediates recruitment of dendritic cells, NK cells, effector Tcells (79) and, in inflamed LNs, subsets of monocytes (80). Both Ox1 and Ox2 also showed upregulated Psme2, a proteosomal component involved in peptide processing for class I antigen presentation (81) (Figure 6C).
While most gene changes were shared, some were preferentially or more dramatically increased in one or the other subset ( Figure 6C). Ccl20, which is selectively expressed by fLEC in resting LNs (Figure 1D), is highly upregulated in Ox1 but not Ox2, maintaining its selective pattern of expression in the SCS floor (Figures 1D, 5B, 6C). Ox1 also displays selective upregulation of Ubiquitin D (Ubd), important in activation of the transcription factor NF-kappa B (82). Oxa controls subset-specific upregulation of monocyte chemoattractants Ccl2 (Ox2) and Ccl5 (Ox1). CCL5 acts through CCR5, which is expressed by multiple adaptive and innate immune subsets; while CCL2 acts through CCR2 in monocyte recruitment (43). Thus, chemotactic regulations in the SCS and medulla diverge. Gene enrichment analysis confirms Ox1 and Ox2 share enrichment in interferon signaling and immune system signatures, but the subsets also display unique ontologies: e.g., Transforming Growth Factor (TGF)-beta regulation of ECM in Ox1 and response to Tumor Necrosis Factor (TNF) in Ox2 ( Figure 6D). Several of the observed changes, including induction of Ccl2, Ccl5, Ccl20, Cxcl9, and Cxcl10 were observed previously in analyses of bulk LN LECs in response to HSV-1 virus (77) or after ovalbumin specific T-cell responses (11). Our results link inflammation-induced transcriptional changes to specific LEC subsets, implying that niche-specific changes coordinate LEC-driven responses in early inflammation.
Interestingly, in both Ox1 and Ox2, significantly more genes are downregulated than upregulated ( Figure 6C). Pathway analysis shows downregulation of genes associated with viral gene expression, ribosomes and mRNA processing ( Figure 6D). This may be a result of the activation of type I interferon signaling (Figure 6D), which can cause widespread downregulation of host and viral transcriptional and translational pathways as a defense against viral replication (83,84). Downregulated genes include Csf1 and Ltb (Figure 6C), which is interesting from the perspective of the known transient loss of sinusoidal LN macrophages in acute inflammation (18,19). Core lymphatic endothelial differentiation genes are also downregulated, including Lyve1 and Cldn5, and subset marker genes are inhibited as well. Ox1 shows reduced expression of fLEC marker genes Cd74 and Bmp2, and Ox2 of Ptx3-LEC marker genes including Stab2 and Tgfbi (Figure 6C). In addition to IFN-induced transcriptional repression, another contributing factor may be transcription factor availability, which can dictate level of gene expression (85): TNF-induced RELA-dependent pathways, which are also induced as part of the Oxa response ( Figure 6D), have been shown to redirect cofactors from superenhancers, repressing cell identity genes in a cell type dependent manner (86).

Common Mouse and Human LN LEC Phenotypes
We have previously reported LEC diversity and described multiple LEC subsets in human LNs (27). Here we compare our mouse LEC samples to a representative sample of human head and neck LN LECs, selected for its quality, high cell number, and inclusion of each of the subsets we identified. We translated human gene names to their mouse homologs, and integrated human and mouse LECs using Canonical Correlation Analysis (CCA) (87), which implements a variant of the mutual nearest neighbor algorithm (MNN) (88). Cross mapping successfully aligned human and mouse LECs, and unsupervised clustering of the aligned cells identified shared but also unique subsets in human LN (Figure 7A). The five subsets common to mouse and human (cLEC, fLEC, Ptx3-LEC, valve-LEC, and Marco-LEC) display high correspondence with human LEC subsets defined previously using unsupervised clustering (27) (LEC I, II, IV, V, and VI, respectively) ( Figure 7A). Specifically, the Marco-LEC subset of medullary EC maps to human LEC VI, which we previously characterized as medullary sinus CD209 + LECs (27). However, mouse Ptx3-LECs map to human LEC IV, which we postulated to be peripheral capillary LECs (27). We discuss this below.

Specialization of the Human SCS Ceiling
Unsupervised clustering aligned the mouse cLECs with one subset of the human cLECs, "(core) cLECs, " which we identified previously as LEC I (27) (Figure 7A), but additional cLEC subsets (cLEC s1, s2, and s3) were unique to the analyzed human LN. Pearson correlation analysis reveals similarity of these subsets with the core cLECs ( Figure 7C). cLEC s3 is most similar to the core cLEC population as indicated by correlation and hierarchical clustering (Figure 7C), UMAP (Figure 7A), and differentially expressed genes (DEGs) (data not shown): we do not discuss this subset further. cLEC s1 corresponds to MFAP4 + SCS LEC (LEC III), which we located overlying the human LN medulla (27). Consistent with this, cLEC s1 clusters close to Ptx3-LECs (Figure 7A), and links other cLECs to Ptx3-LECs in trajectory space (see below Figure 10A). cLEC s2, which we did not segregate previously, uniquely expresses high levels of Hairy/enhancer-of-split related with YRPW motif protein (HEY1), the chemokine CCL2 (Figure 7C), and E-Selectin (SELE) (data not shown), suggestive of an activated state. Both E-selectin and CCL2 are induced by inflammation in endothelial cells and regulate leukocyte recruitment (89). Based on HEY1 and SELE expression, the subset is identifiable in 3 of the 6 human LN samples we studied previously (27) (data not shown). The findings suggest greater heterogeneity and specialization in SCS ceiling LEC in humans than in the resting specific-pathogenfree (SPF) mice studied here. The subset specialization of human cLEC may relate to local differences in immune environments, or to the more complex architecture of the human LN. Unlike the mouse, human cLECs participate in invaginations of the capsule, known as trabecular sinuses (4), which may experience more turbulent flow of the incoming lymph and hence variation in shear stress.

Comparisons of Human and Mouse LEC Differentially Expressed Genes (DEG)
To evaluate similarities between species, we compared DEGs of mouse and human LEC subsets using a gene overlap score, defined as the ratio of the number of shared DEGs to the number of genes expected to be shared based on random chance, for each combination of mouse and human subsets. In all instances, the highest overlap scores are seen between corresponding subsets ( Figure 7B). Based on overlap scores, cLECs are more conserved than the more immunologically active subsets fLECs, Marco-LECs, and Ptx3-LECs. The floor and medullary sinus subsets showed less conserved DEG profiles, likely reflecting evolutionary pressure in response to pathogens, contrasting with conservation of structural functions of cLECs. Gene set enrichment analysis based on conserved genes has the potential to highlight core functions of the LEC subsets ( Figure 7D). Shared cLEC genes are involved in focal adhesion, fluid shear stress response and ECM interaction, consistent with their structural role and association to the capsule. Ptx3-LECs are enriched in genes for endocytosis and lysosome, which could relate to the high expression of scavenger receptors (e.g., Lyve1, Stab1, Stab2). Marco-LECs are enriched for lysosomal genes, and for genes of the JAK-STAT signaling pathway as well as the coagulation and complement cascades. fLECs display high enrichment for inflammatory pathways, including JAK-STAT and Nuclear Factor (NF)-kappa B, and are enriched for genes involved in antigen processing and presentation, supporting their immunomodulatory properties. A number of mouse key subset specific marker DEGs are shared in the human, including fLECs: Bmp2, Ccl20, Cd74, and Csf1; cLECs: Ackr4, Bmp4, and Cav1; Marco-LECs: Marco, Clec4g, and C2; and Ptx3-LECs: Ptx3, Nrp2, and Flt4 (Vegfr3). We illustrate expression of these and other select shared marker genes in Figure 8.
As noted above, mouse Marco-LECs align with LEC VI, which we identified previously as the major medullary sinus subset in human LNs (27). Unexpectedly, mouse Ptx3-LECs align with a subset (LEC IV) that we previously related to capillary LECs based on their enrichment for expression of Ccl21 and Lyve1 (27), gene markers of peripheral capillary lymphatic vessels (57,58). However, these genes are also expressed by mouse Ptx3-LECs (Figures 1D, 5B, 8), and as noted earlier, their expression likely reflects the parallels with capillary LEC in morphology (blind ends, loose EC junctions) and function (recruiting fluid and lymphocytes into lymph) (4). Supporting a medullary identity, Ptx3-LECs and LEC IV share high levels the sinusoidal endothelial marker Stab2/STAB2 (46), shown to be expressed by medullary sinuses in human LNs (90). Importantly, human Ptx3-LECs (LEC IV) also share expression of PTX3 and lack PD-L1 (CD274) expression, similar to their mouse counterpart; and they express the inter-alpha-trypsin inhibitor gene family member ITIH3, functionally related to the mouse Ptx3-LEC marker Itih5 (64) (Figure 8).
Ptx3-LECs in both humans and mice are distinguished from Marco-LECs by higher expression of the glycoprotein and scavenger receptor CD36, also known as Fatty Acid Translocase (FAT) (Figure 9A). Staining of CD36 in human head and neck LNs identified capillary-like, CD36 high LYVE-1 + lymphatic cords, negative for MARCO and the Marco-LEC [LEC VI (27)] marker CLEC4M (Figures 9B,C). These CD36 high lymphatic sinuses were found either as isolated cords in the paracortex (Figures 9B,C) or as extended sprouts from MARCO + CLEC4M + medullary sinuses (Figure 9D), similar to mouse Ptx3-LECs which connects to Marco-LECs in perifollicular regions (Figure 2). Human Ptx3-LECs are also distinguished from Marco-LECs by high expression of PDPN and CCL21, as predicted by the transcriptional profiles (Figure S8), while both medullary subsets express LYVE-1 and STAB2 (Figure S8). Thus, cross-species comparison of single-cell profiles (Figure 7) and in situ analysis (Figure 9) support that human and mouse share two distinct niches of paracortical and medullary sinus LECs: MARCO + LECs, which correspond to the previously published CD209 + medullary sinus LEC subset (LEC VI) (27), and Ptx3-LECs, PDPN high , CD36 high paracortical and medullary sinus LECs, corresponding to LEC IV (27) in our earlier classification of human LN LECs. Notably, human and mouse LN LECs also share the in situ physical transition between Marco-LECs and Ptx3-LECs, as predicted in trajectory analysis (Figures 1C, 10A, discussed below).

Species-Specific Gene Expression
A number of genes with homologs in mouse and human are not conserved in expression, or display different patterns of subset selectivity. Here we highlight select examples for discussion ( Figure S9). As noted earlier, the mouse fLEC markers Madcam1 and Msr1 are poorly expressed by human LECs (27) (data not shown). ACKR1, also known as DARC (Duffy Antigen Receptor for Chemokines), marks human fLECs and Marco-LECs, but Ackr1 has very low expression in mouse LN LECs, without clear subset selectivity (Figure S9B). Consistent with this, endothelial ACKR1 in mouse is restricted to venular blood vessels, with only sparse and low detection in LECs (91). ACKR1 is a chemokine "interceptor, " which can serve as a sink for a large range of inflammatory chemokines (92). Its expression could reflect a greater need to moderate inflammatory chemokines in the human. Alternatively, it may facilitate transport of chemokines across the lymphatic endothelium in the human LN, as ACKR1 has been shown to shuttle chemokines across endothelial cell layers (93). Several ACKR1 ligands are expressed by human fLECs including CXCL3 and CXCL5 (27).
Human but not mouse LN LECs also display high expression of IL6 in fLECs (Figure S9B), likely reflecting a higher inflammatory basal state in human, especially compared to our SPF mice. Human Ptx3-LECs express MMP2 and LOX, genes missing in mouse LECs or expressed in different subsets (i.e., Lox expressed weakly in cLECs in mouse) ( Figure S9B). Both Matrix metalloproteinase 2 (MMP2) and Lysyl oxidase (LOX) can contribute to ECM modulation (94,95), suggesting that, although matrix interplay is conserved, the specific mechanisms of matrix interaction in this population of LECs may have diverged. Mouse-specific LN LEC DEGs include Apolipoprotein E (ApoE) in cLECs and medullary sinus populations and Regulator of G-protein signaling 4 (Rgs4) in fLECs ( Figure S9B). Carboxypeptidase E (Cpe) and Carbonic anhydrase 8 (Car8) are examples of genes with different expression pattern across mouse and mouse LN LEC subsets ( Figure S9B). Since scRNA-seq is often unable to detect low abundance transcripts, the apparent lack of expression of a gene must be interpreted with caution: the expression pattern of the genes mentioned here have been observed in each of our samples.
In addition, CD209 and CLEC4M, which lack orthologs in mouse (27), are specific for human Marco-LECs. We described them previously as medullary sinus markers (27). Human IL32 and mouse Glycam1 also lack orthologs. IL-32 potently amplifies inflammatory responses by induction of multiple cytokines (96); it is highly expressed by human Ptx3-LECs (data not shown). GLYCAM1, a secreted mucin that on high endothelial venules is decorated with glycotopes for leukocyte selectins (97), is selectively expressed by mouse fLEC (Figure 5B and Figure S4).

The Continuum of Human LEC Phenotypes in Trajectory Space
We have focused to this point on comparing "subsets" of human and mouse LECs, but as illustrated for the mouse, LECs exist in a phenotypic continuum that may reflect physical alignments (spatial transitions), developmental sequences, or both. To gain further insight into LEC diversity and zonation within human LNs, we ran tSpace on combined LN LECs from our six previously published human samples (27). In the tSpace projections in Figure 10A, cell identities are determined based on correlation with the core subsets in our index human LN (left) or based on unsupervised clustering as previously (27) (right). Alignment of human LECs in trajectory space reveals both similarities and differences to the mouse. Shared subsets are aligned in the same order as in the mouse, with links from Valve-LECs to cLECs, Ptx3-LECs, Marco-LECs and finally fLECs (Figure 10A) (compare with mouse alignments in Figures 1B,C). However, in the human, the cellular "bridge" between cLECs and fLECs is highly populated with cells. The human-specific cLEC subsets also show interesting alignments. cLEC s2 appears closely linked to the bridge. cLEC s1, which we previously identified as SCS ceiling LEC overlying the medulla (LEC III) (27), extends toward and links prominently to Ptx3-LECs (see arrows): Cells from both "subsets" within this Ptx3-LECs-to-cLEC s1 trajectory express PTX3 itself, and are highly enriched in CCL21 (Figure 10B), suggesting a relationship of PTX3 + medullary sinuses and the perihilar ceiling. Supporting this notion, RNAscope of a neonatal human LN showed extensive PTX3 in the medulla of the hilar region (data not shown). Co-expression of CCL21 and SPHK1 (enzyme required for S1P production) support a role of Ptx3-LECs in lymphocyte egress. Heterogeneity in ANGPT2 also exists along the Ptx3-LEC paths ( Figure 10B).
The alignments outlined here are not an artifact of the high cell number or integration of LEC from different LNs, because the same patterns and linkages are seen in independent tSpace projections of our index LN sample (not shown). Human LN LECs thus display numerous intermediate phenotypes and complex zonation that may reflect the complexity of human LN architecture observed histologically.

Summary and Outlook
The LN lymphatic vasculature provides a complex lymphatic vascular bed, adapted to a unique microenvironment where cooperation between stromal cells and immune cells forms the basis for effective immune cell interaction and activation.
Endothelial cells in this environment not only need to provide a framework for the structural organization of the organ, but also need to be able to adapt to constant changes induced by immunological stimuli and organ expansion in LN hypertrophy, and to support tolerogenic immune reactions in homeostasis.
We define five LN LEC subsets in mouse: valve LECs, the SCS Ackr4 + cLECs, the immune active Ccl20 + fLECs, as well as Ptx3-LECs and Marco-LECs, two paracortical and medullary sinus subsets that were not recognized previously. Single-cell gene profiles indicate niche-specific functional specialization of these subsets with distinct pathways for pathogen interactions and matrix modeling. Interestingly, cross-species mapping with human LN LECs shows that both subsets as well as their respective functions are conserved, which allows us to redefine the subset previously identified as candidate capillary LECs (LEC IV) (27) in human LN to paracortical sinus Ptx3-LECs. PTX3 and MARCO, gene markers chosen here to designate the two distinct medullary LEC subsets, are conserved both in their expression patterns among LEC subsets and also in protein sequence between humans and mice (92% homology PTX3 and 69% MARCO), suggesting that they may have conserved functions as well.
Recently developed algorithms have made important advances in solving the complex problem of integrating different datasets (87)(88)(89)(90)(91)(92)(93)(94)(95)(96)(97)(98)(99). As shown here, they not only perform well in combining sample replicates by removing "batch effects, " but can also identify and map similar subsets of cells across species barriers. Mutual nearest neighbors and CCA algorithms not only mapped human LEC subsets to mouse counterparts, but also uncovered species-specific LN LEC subset compositions, revealing a greater diversity of SCS ceiling LEC in human. We can look forward to continuing advances in computational approaches to integrating and mining scRNA-seq data. Particularly exciting is the power of trajectory inference to recapitulate or predict the organization of endothelial cells within the complex vascular networks. While trajectory analysis has been shown to model sequences of cell phenotypes in development ("pseudotime") (30), our results both in lymphatics (here) and blood vascular EC studies (31) show that computed alignments of cells in trajectory space reflect the tissue architecture (spatial organization) and physical relationships between endothelial cells in situ with surprising faithfulness. This correspondence implies, as an approximation, that endothelial cell phenotypes progress in a gradual and orderly fashion within the linear arrangements (as in vascular tubes) or sheets (as in the SCS and sinuslining LECs in the LN) that make up the endothelium. In essence, trajectories provide a computational roadmap for mapping gene expression to the vascular endothelium. The results underscore the diversity of endothelial cells as a continuum, punctuated by concentrations of particular phenotypes or niches that are identified as subsets. Whether the progression of phenotypes and zonation among LN LECs reflects malleable LEC responses to local niche factors, retention of a programmed developmental response, or both, remains to be determined.
Our studies here demonstrate the power of single-cell profiling to illuminate the biology of the vascular endothelium, and the promise it holds to revolutionize our understanding of conserved and species-specific regulation of the vasculature and its responses in physiology and human disease.

Mice
Male and female 6-8-week-old BALB/cJ or 8-12-week-old C57BL/6J peripheral LNs (inguinal, axillary and brachial) (processed by 10x Genomics workflow), or 20-week-old Prox1-GFP/C57BL/6J (33) female inguinal LNs (processed by SMART-seq2 workflow) were used for scRNA-seq. Mice were bred and maintained in the animal facilities at Veterans Affairs Palo Alto Health Care System, accredited by the Association for Assessment and Accreditation of Laboratory Animal. In addition, mice were held at a Specific Pathogen Free (SPF) facility Uppsala University, and experimental procedures were approved by the local animal ethics committee at the Uppsala County Court (6009/17).
Count data were processed with the Seurat package (v3.1.0) (87, 99). For quality control, genes that were expressed in fewer than three cells and cells that expressed fewer than 100 genes were excluded from analysis. Raw counts were log normalized, and 2000 most variable genes were identified based on a variance stabilizing transformation. To correct for batch differences, variable gene sets were used to align multiple datasets for joint analysis, using the Canonical Correlation Analysis (CCA) method within the "FindIntegrationAnchors" and "IntegrateData" functions of the Seurat package, which implements a variant of the mutual nearest neighbors (MNN) algorithm. Correspondences between single cells across datasets were identified and multiple datasets were transformed into a shared space. Principal Component Analysis (PCA) dimensionality reduction was performed using the variable gene sets. Cell clusters were determined using a Shared Nearest Neighbor (SNN) modularity optimizationbased clustering algorithm of the Seurat "FindClusters" function, and were visualized with t-distributed Stochastic Neighbor Embedding (tSNE) (102) or Uniform Manifold Approximation and Projection (UMAP) (103). Contaminating pericyte, immune and blood endothelial cells were removed by supervised gating on the tSNE plot. To recover gene-gene relationships that are lost due to dropouts, we imputed missing gene expression data from log normalized count data using an in-house customization (https://github.com/kbrulois/magicBatch) of the MAGIC (Markov Affinity-based Graph Imputation of Cells) algorithm with optimized parameters (t = 2, k = 9, ka = 3) (104). Imputed data were used for visualization of single-cell gene expression in violin plots and heatmaps (105), as well as for trajectory analyses.

Nearest Neighbor Alignments in Trajectory Space
tSpace was used to model the nearest neighbor relationships as well as transitional zones between LEC subsets (30). Mapping of cells in trajectory space allows unsupervised reconstruction and exploration of complex developmental sequences, without implementation of cell clustering or manifold distortion for graphical expedience. Batch effects were removed using the "fastMNN" function of the batchelor package (v1.0.1) (88) and variable genes were identified using the scran package (106). For mouse LEC, we used tSpace with default parameters (T = 100, K = 20, 5 graphs, 20 way points, Euclidean distance metric) and imputed expression values of the top 800 variable genes (Figures 1C, 6A, right panel) or batch-corrected low-dimensional coordinates (top 50 coordinates) (Figure 6A, left panel). The trajectory matrices were visualized in low dimensional space using PCA (Figures 1C, 6A, right panel) or UMAP ( Figure 6A, left panel) within the tSpace package. For combined human LEC samples, imputed data of LEC from 6 LNs were batch-corrected as above. Twenty PCs from the 1000 variable genes were calculated, the loadings were adjusted to a minimum value of 0 by addition of a constant and used as input to the tSpace algorithm (T = 200, K = 25, 5 graphs, 10-way points, cosine distance metric).

Differential Gene Expression and Gene Enrichment Analysis
Differential gene expression analysis was performed using the negative binomial generalized linear model within the "FindMarkers" function of Seurat on log normalized count data (p < 0.01, fold change > 1.2). To identify subset-defining genes, single-cell transcriptomes of cells in the subset of interest were compared to those of all other cells within the same sample. Only upregulated DEGs were considered, unless otherwise specified. Gene enrichment analysis of DEGs with Gene Ontology (107,108), Kyoto Encyclopedia of Genes and Genomes (109,110), and BioPlanet (111) databases was performed using Enrichr (112,113).

Correlation Analysis
To determine similarities between LEC subsets, mean gene expression values were calculated for each subset. The top 1000 most variable genes across all subsets were used for pairwise Pearson correlation analysis. Hierarchical clustering was performed using the average linkage method.

Cross Species Single-Cell Transcriptome Analysis
Human HNLN1 dataset (GSE124494) (27) was used for the integrated LEC profiling. Human gene names were converted to their mouse homologs using the biomaRt package (v2.40.4) (114). scRNA-seq datasets were integrated in Seurat for unsupervised clustering. The top 100 most upregulated DEGs for each subset, as ranked by the fold change of gene expression in the subset relative to other LEC combined, were determined for human LEC and separately for mouse. An overlap score was defined as the number of DEGs common to one human subset and one mouse subset divided by the number of genes expected to be shared by the two subsets by chance. For each pair of mouse subset and human subsets, overlap score = n / (A × B / N), where n = number of observed overlapping genes between the top 100 DEGs of the human subset and the mouse subset, A = number of DEGs considered in the mouse subset (100), B = number of DEGs considered in the human subset with mouse homologs, and N = total number of genes detected in both mouse and human LEC, 13,458.

Reanalysis of Human LEC Datasets
For the combined analysis of human LEC samples (27), imputed gene expression data were batch-corrected and used for trajectory analysis as above. LEC in each of the additional five human samples were classified by correspondence to the index human subsets that mapped with mouse subsets as follows: reference subset mean gene expression was generated from the index HNLN1 human dataset, using core cells of each major subset that cross-mapped with mouse subsets, as well as manually gated bridge cells that link fLEC and cLEC in tSpace projections. Cells in the other human samples were classified by Pearson correlation using the 1000 most variable genes in the reference set.

Immunostaining of Mouse LNs
Inguinal and popliteal mouse LNs were harvested from Prox1-GFP mice (33) and fixed in 0.8% paraformaldehyde (PFA) for 12 h at 4 • C. After fixation, the LNs were placed in sucrose: 25% for 2 days, 50% for 2 days before embedding in OCT media (HistoLab), then snap frozen on dry ice. Frozen tissues were cryo-sectioned at a thickness of 8 µm and stored at −80 • C. For fresh frozen tissue, LNs were harvested from wild-type C57BL/6 mice and cleaned in ice cold PBS, embedded in OCT media, and snap frozen on dry ice. For immunostaining, the sections were hydrated in Phosphate-buffered saline (PBS) and blocked with 10% donkey serum (Sigma) diluted in PBS for 20 min. After blocking, the sections were incubated with primary antibodies diluted in blocking buffer overnight at 4 • C. Thereafter the sections were washed in PBS with 0.1% TritonX100 (Sigma) (PBSTX) and incubated with secondary antibodies diluted in PBSTX for 1 h at RT. The sections were counterstained with 4 ′ ,6diamidino-2-phenylindole (DAPI) following additional washing in PBSTX and mounted in ProLong Gold Antifade Mountant (Thermofisher Scientific).

Immunostaining of Human LNs
Human head and neck tumor-free LNs from cancer patients were received from the hospital immediately after the surgery, and embedded in OCT compound (Sigma) and frozen on dry ice. The collection was done under the license ETMK: 132/2016. A written informed consent was obtained from each individual donating tissue. The samples were kept anonymous and used with the permission of the Ethical Committee of Turku University Hospital. The LNs were sectioned at a thickness of 6 µm with a cryostat, and fixed with acetone at −20 • C. The sections were incubated with 10% FCS for blocking, incubated with the primary antibodies diluted in 0.5% BSA in PBS overnight at 4 • C. Thereafter, they were incubated with the secondary antibodies for 2 h at room temperature. Sections were washed with PBS and mounted with ProLong Gold Antifade Mounting medium with DAPI (Thermofisher Scientific).

RNA in situ Hybridization
Inguinal mouse LNs were harvested from wild-type C57BL/6 mice following fixation with 2% PFA by heart perfusion. The tissues were embedded in OCT and snap frozen on dry ice. The tissues were cryosectioned at a thickness of 14 µm, dried at −20 • C for 1 h and stored at −80 • C. In situ hybridization (ISH) was performed using RNAscope Multiplex Fluorescent kit according to the manufacturer's instructions (Advanced Cell Diagnostics). Briefly, the sections were fixed in ice-cold 4% PFA for 15 min, rinsed in PBS and dehydrated with increasing concentrations of ethanol: 50, 70% and absolute ethanol for 5 min each. The sections were dried at RT and treated with protease IV for 15 min and rinsed in PBS. Thereafter the sections were incubated with the mouse probes: Claudin5-C3, Stabilin1-C1, Bmp2-C2, Bmp4-C2, Ptx3-C1, Prox1-C2 in ACD HybEZ II hybridization system (Advanced Cell Diagnostics) at 40 • C for 2 h. The remainder of the assay protocol was implemented following manufacturer's statement. For detection of the combinations Ptx3-Prox1 and Bmp4-Claudin5 we performed the RNAscope assay with amplification V2 Kit (Advanced Cell Diagnosis). The sections were counterstained with DAPI.

Imaging
Images were obtained with Vectra Polaris TM Automated Quantitative Pathology Imaging System (Akoya Biosciences) and LSM 700 or LSM 880 confocal microscopes (Zeiss). Confocal objectives: Air objective plan apo N/A 0.45 (10x magnification); air objective plan apo N/A 0.8 (20x magnification); air objective plan apo N/A 0.95 (40x magnification). Analyses were performed with ImageJ software.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethical Committee of Turku University Hospital. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by the Animal Ethics Committees of Veterans Affairs Palo Alto Health Care System and Uppsala County Court.