Sphingosine 1-phosphate-regulated transcriptomes in heterogenous arterial and lymphatic endothelium of the aorta

Despite the medical importance of G protein-coupled receptors (GPCRs), in vivo cellular heterogeneity of GPCR signaling and downstream transcriptional responses are not understood. We report the comprehensive characterization of transcriptomes (bulk and single-cell) and chromatin domains regulated by sphingosine 1-phosphate receptor-1 (S1PR1) in adult mouse aortic endothelial cells. First, S1PR1 regulates NFκB and nuclear glucocorticoid receptor pathways to suppress inflammation-related mRNAs. Second, S1PR1 signaling in the heterogenous endothelial cell (EC) subtypes occurs at spatially-distinct areas of the aorta. For example, a transcriptomically distinct arterial EC population at vascular branch points (aEC1) exhibits ligand-independent S1PR1/ß-arrestin coupling. In contrast, circulatory S1P-dependent S1PR1/ß-arrestin coupling was observed in non-branch point aEC2 cells that exhibit an inflammatory gene expression signature. Moreover, S1P/S1PR1 signaling regulates the expression of lymphangiogenic and inflammation-related transcripts in an adventitial lymphatic EC (LEC) population in a ligand-dependent manner. These insights add resolution to existing concepts of endothelial heterogeneity, GPCR signaling and S1P biology.


Introduction
Sphingosine 1-phosphate (S1P), a circulating lipid mediator, acts on G protein-coupled S1P receptors (S1PRs) to regulate a variety of organ systems. S1PR1, abundantly expressed by vascular endothelial cells (ECs), responds to both circulating and locally-produced S1P to regulate vascular development, endothelial barrier function, vasodilatation and inflammation . S1P binding to S1PR1 activates heterotrimeric G ai/o proteins, which regulate downstream signaling molecules such as protein kinases, the small GTPase RAC1, and other effector molecules to influence cell behaviors such as shape, migration, adhesion and cell-cell interactions. Even though S1PR1 signaling is thought to evoke transcriptional responses that couple rapid signal transduction events to long-term changes in cell behavior, such mechanisms are poorly understood, especially in the vascular system. Subsequent to activation of G ai/o proteins and RAC1, the S1PR1 C-terminal tail gets phosphorylated and binds to ß-arrestin, leading to receptor desensitization and endocytosis (Liu et al., 1999;Oo et al., 2007). While S1PR1 can be recycled back to the cell surface for subsequent signaling, sustained receptor internalization brought about by supra-physiological S1P stimulation or functional antagonists that are in therapeutic use leads to recruitment of WWP2 ubiquitin ligase and lysosomal/ proteasomal degradation of the receptor (Oo et al., 2011). Thus, ß-arrestin coupling down-regulates S1PR1 signals. However, studies of other GPCRs suggest that ß-arrestin coupling can lead to biased signaling distinct from G ai/o -regulated events (Wisler et al., 2018). Distinct transcriptional changes brought about by G ai/o -and ß-arrestin-dependent pathways are not known.
Our understanding of GPCR signaling in vivo, particularly with respect to direct transcriptional targets and spatial specificity of signaling, is limited. To address this, Kono et al. (2014) developed S1PR1 reporter mice (S1PR1-GS mice) which record receptor activation at single-cell resolution (Kono et al., 2014;Barnea et al., 2008). S1PR1-GS mice harbor one wild-type S1pr1 allele and one targeted knock-in allele, which encodes S1PR1-tTA and ß-arrestin-TEV protease fusion proteins separated by an internal ribosome entry sequence (Kono et al., 2014). Breeding the S1pr1 knock-in allele with the tTA-responsive H2B-GFP allele generates an S1PR1-GS mouse. In S1PR1-GS mice, the b-arrestin-TEV fusion protein triggers release of tTA from the C terminus of modified S1PR1 when barrestin-TEV and S1PR1-tTA are in close proximity. Free tTA enters the nucleus and activates H2B-GFP reporter gene expression. Since S1PR1-GS mouse embryonic fibroblast cells respond to S1P with an EC 50 of 43 nM (Kono et al., 2014;Lee et al., 1998), the S1PR1 reporter system accurately reports S1PR1 activation. Indeed, structural and functional analyses of other GPCRs suggest that the C-terminal phosphorylation patterns determine the strength of ß-arrestin binding, which was accurately predicted by the ß-arrestin coupling (Zhou et al., 2017). The in vivo half-life of H2B-GFP protein is~24 days in hair follicle stem cells (Waghmare et al., 2008). Therefore, GFP expression in this reporter mouse represents the cumulative record of S1PR1 activation in vivo.
We previously showed that high levels of endothelial GFP expression (i.e. S1PR1/ß-arrestin coupling) are prominent at the lesser curvature of the aortic arch and the orifices of intercostal branch points (Galvani et al., 2015). In addition, inflammatory stimuli (e.g. lipopolysaccharide) induced rapid coupling of S1PR1 to ß-arrestin and GFP expression in endothelium in an S1P-dependent manner (Kono et al., 2014). These data suggest that the S1PR1-GS mouse is a valid model to study GPCR activation in vascular ECs in vivo.
To gain insights into the molecular mechanisms of S1PR1 regulation of endothelial transcription and the heterogenous nature of S1PR1 signaling in vivo, we performed bulk transcriptome and open chromatin profiling of GFP high and GFP low aortic ECs from S1PR1-GS mice. We also performed transcriptome and open chromatin profiling of aortic ECs in which S1pr1 was genetically ablated (S1pr1 ECKO) (Galvani et al., 2015). In addition, we conducted single-cell (sc) RNA-seq of GFP low and GFP high aortic ECs. Our results show that S1PR1 suppresses the expression of inflammation-related mRNAs by inhibiting the NFkB pathway. Second, the high S1PR1 signaling ECs (GFP high cells) are more similar to S1pr1 ECKO ECs at the level of the transcriptome. Third, scRNA-seq revealed eight distinct aorta-associated EC populations including six arterial EC subtypes, adventitial lymphatic ECs, and venous ECs, the latter likely from the vasa vasorum. S1PR1 signaling was highly heterogenous within these EC subtypes but was most frequent in adventitial LECs and two arterial EC populations. Immunohistochemical analyses revealed spatio-temporal regulation of aortic EC heterogeneity. In lymphatic ECs of the aorta, S1PR1 signaling restrains inflammatory and immunerelated transcripts. These studies provide a comprehensive resource of transcriptional signatures in aortic ECs, which will be useful to further investigate the multiple roles of S1P in vascular physiology and disease.

Results
Profiling the transcriptome of GFP high and GFP low mouse aortic endothelium To examine S1PR1/ß-arrestin coupling in the aorta, we used the previously described S1PR1-GS (Kono et al., 2014) mouse strain. Mice heterozygous for the knock-in allele (S1pr1 ki/+ ) are born at the expected Mendelian frequency (Figure 1-figure supplement 1A) and do not show phenotypic abnormalities (Figure 1-figure supplement 1B and C). However, homozygous mice (S1pr1 ki/ki ) showed an~2 fold reduction in circulating lymphocytes and~2 fold increase in lung vascular leakage of Evans Blue dye suggesting that hypomorphism of the fusion S1pr1 allele in the signaling mouse. Therefore, all subsequent experiments were performed using heterozygous S1pr1 ki/+ mice harboring one allele of the H2B-GFP (Tumbar et al., 2004) reporter gene, which do not exhibit S1pr1 hypomorphic phenotypes. S1PR1 expression in aortic endothelium is relatively uniform (Galvani et al., 2015). However, S1PR1 coupling to ß-arrestin, as reported by H2B-GFP expression in S1PR1-GS mice, exhibits clear differences in specific areas of the aorta. For example, thoracic aortae of S1PR1-GS mice show high levels of GFP expression in ECs at intercostal branch points (Galvani et al., 2015) but not in ECs of control (S1pr1 +/+ ) mice harboring only the H2B-GFP reporter allele ( Figure 1A), confirming that GFP expression in aortic ECs is dependent on the S1pr1 knock-in allele. The first 2-3 rows of cells around the circumference of branch point orifices exhibit the greatest GFP expression ( Figure 1A). In addition, heterogeneously dispersed non-branch point GFP+ ECs were also observed, including at the lesser curvature of the aortic arch ( Figure 1A). Areas of the aorta that are distal (>~10 cells) from branch points, as well as the greater curvature, exhibit relatively low frequencies of GFP+ ECs ( Figure 1A). GFP+ mouse aortic ECs (MAECs) are not co-localized with Ki-67, a marker of proliferation, suggesting that these cells are not actively cycling ( Figure 1A). However, fibrinogen staining was frequently co-localized with GFP+ MAECs, suggesting that ß-arrestin recruitment to S1PR1 was associated with increased vascular leak ( Figure 1A). These findings suggest sharp spatial differences in S1PR1 signaling throughout the normal mouse aortic endothelium.
Differential expression analysis identified 1,103 GFP high -enriched and 1,042 GFP low -enriched transcripts (p-value<0.05) ( Figure 1C and Figure 1-figure supplement 2D; see also Supplementary file 1). In contrast, S1pr1 ECKO MAECs showed fewer differentially expressed genes (DEGs), with 258 up-and 107 down-regulated transcripts ( Figure 1C and Figure 1-figure supplement 2E; see also Supplementary file 1). Intersection of these two sets of DEGs showed that only 9.5% (204 transcripts) were common ( Figure 1C and Supplementary file 1), suggesting that the majority (~90%) of transcripts that are differentially expressed in MAECs from S1PR1-GS mice are not regulated by S1PR1 signaling. Rather, S1PR1/ß-arrestin coupling correlates with heterogenous EC subtypes in the mouse aorta.
Among the 204 common DEGs, 151 were both S1pr1 ECKO up-regulated and enriched in the GFP high population ( Figure 1C). In contrast, much lower numbers of transcripts were found in the intersection of GFP high and S1pr1 ECKO down-regulated (seven transcripts), GFP low and S1pr1 ECKO up-regulated (eight transcripts) and GFP low and S1pr1 ECKO down-regulated (38 transcripts) Figure 1. High S1PR1/ß-arrestin coupling in normal mouse aortic endothelium exhibits transcriptomic concordance with S1PR1 loss-of-function. (A) H2B-GFP control and S1PR1-GS mouse thoracic aorta whole-mount en face preparations. Representative images from different regions of the aorta are presented and show H2B-GFP (GFP), VE-Cadherin (VEC) or CD31, Ki67 (N = 3) or Fibrinogen (N = 2) immunostaining. Scale bars are 50 mM. (B) FACS gating scheme used for isolation of GFP high and GFP low MAECs showing the uncompensated CD31-PE and GFP channels. S1pr1 ECKO and WT MAECs were isolated using the GFP low gate of this scheme (see  Figure 1C) (detailed gene set overlaps are provided in Supplementary file 1). We computed the statistical significance of these gene set overlaps using the GeneOverlap R package (Shen et al., 2013). The GFP high :S1pr1 ECKO-up overlap was significant (151 transcripts, p-value=3.80E-126), as was the GFP low :S1pr1 ECKO-down overlap (38 transcripts, p-value=1.3E-22), and the other two tested overlaps were not significant (p-value>0.3) (Figure 1-figure supplement 2F). These data suggest that the transcriptome of MAECs exhibiting high S1PR1/ß-arrestin coupling (GFP high ) is more similar to that of S1pr1 ECKO.
We used Ingenuity Pathway Analysis (IPA, Qiagen) to examine biological processes regulated by S1PR1 signaling and loss of function in MAECs. Transcripts involved in inflammatory processes were prominently up-regulated in both GFP high and S1pr1 ECKO MAECs ( Figure 1D; see also Supplementary file 2). For example, positive tumor necrosis factor (TNF)-a, lipopolysaccharide, and interferon-g signaling were observed in both S1pr1 ECKO and GFP high MAECs. In contrast, a negative glucocorticoid signature was observed in these cells ( Figure 1D). Examples of differentially regulated transcripts are chemokines (Ccl2, Clc5, Ccl7, Ccl21c), cytokines (Il33, Il7), inflammatory modulators (Irf8, Nfkbie, Tnfaip8l1) and cyclooxygenase-2 (Ptgs2) ( Figure 1E and Figure 1-figure supplement 2D and E). This suggests that S1PR1 suppresses inflammatory gene expression in mouse aortic endothelium. We noted that transcripts in the TGFß signaling pathway (Thbs1, Smad3, Bmpr1a, Col4a4, Pcolce2) were prominently down-regulated in the GFP high population ( Figure 1D and Figure 1-figure supplement 2D). Furthermore, both GFP high and S1pr1 ECKO MAECs were enriched with Lyve1, Flt4, and Ccl21c transcripts, which encode proteins with well-defined roles in lymphatic EC (LEC) differentiation and function (Ulvmar and Mäkinen, 2016; Figure 1E, Figure 1figure supplement 2G). Taken together, these data suggest that S1PR1 represses expression of inflammatory genes in aortic endothelium and that GFP high MAECs include aorta-associated LECs and are heterogeneous.

Chromatin accessibility landscape of MAECs
We used the assay for transposase-accessible chromatin with sequencing (ATAC-seq) (Buenrostro et al., 2013) to identify putative cis-elements that regulate differential gene expression between GFP high versus GFP low and S1pr1 WT versus S1pr1 ECKO MAECs. ATAC-seq utilizes a hyper-active Tn5 transposase (Adey et al., 2010) that simultaneously cuts DNA and ligates adapters into sterically unhindered chromatin. This allows for amplification and sequencing of open chromatin regions containing transcriptional regulatory domains such as promoters and enhancers. After alignment, reads from three experiments were trimmed to 10 bp, centered on Tn5 cut sites, then merged. These merged reads were used as inputs to generate two peak sets (MACS2, FDR < 0.00001) of 73,492 for GFP low MAECs and 65,694 for GFP high MAECs ( Figure 2A). MAECs isolated from WT and S1pr1 ECKO mice harbored 93,859 and 76,082 peaks, respectively ( Figure 2A). Peaks were enriched in promoter and intragenic regions (Figure 2-figure supplement 1A). We noted that the Cdh5 gene exhibited numerous open chromatin peaks, while Gata1 was inaccessible ( Figure 2-figure supplement 1B). Furthermore, we observed a global correlation between chromatin accessibility and mRNA expression for all 20,626 annotated coding sequences (CDS's) in the NCBI RefSeq database (Figure 2-figure supplement 2). These data suggest that our ATAC-seq data is of sufficient quality for detailed interrogation. of GFP high vs. GFP low (top) and S1pr1 ECKO vs. WT (bottom) MAEC comparisons. Activation Z-scores and P-values are indicated for each selected factor (see Supplementary file 2). (E) Expression heatmaps (row Z-scores) of the 159 S1pr1 ECKO up-regulated transcripts also differentially expressed between GFP high and GFP low MAECs. Values represent individual replicates from comparison of S1pr1 ECKO vs WT (left) and GFP high vs GFP low (right) and selected transcripts are labeled. For both RNA-seq experiments, three cohorts of mice were used for MAEC isolation and downstream statistical analyses (N = 3 for each experiment). The online version of this article includes the following figure supplement(s) for figure 1: Figure supplement 1. S1PR1-GS mice with a single S1pr1-knockin allele are phenotypically normal. Figure supplement 2. RNA-seq quality control and differential gene expression between GFP high and GFP low and S1pr1 ECKO and WT MAECs. Figure 2. S1PR1 loss-of-function and high levels of ß-arrestin coupling are associated with an NFkB signature in open chromatin. (A) Venn diagrams illustrating all peaks (FDR < 0.00001) identified after analysis of individual ATAC-seq replicates of GFP low , GFP high , S1pr1 ECKO, and WT MAECs, and subsequent merging of these peaks into a single consensus peak set. (B) Volcano plots of all peaks for both the GFP high vs GFP low and S1pr1 ECKO vs WT MAECs. Three individual experiments were performed for GFP low vs GFP high and S1pr1 ECKO vs WT comparisons (N = 3). Differentially accessible peaks (DAPs) were determined using edgeR (FDR < 0.05, see Materials and methods) (colored dots). (C-F) Transcription factor (TF) binding motif enrichment analysis of DAPs. The DAPs were input to the HOMER 'findMotifsGenome.pl' script (see also Supplementary file 2) and observed vs expected frequencies of motif occurrances were plotted. (G-H) Graphs showing ATAC signal at predicted TF Figure 2 continued on next page Differential chromatin accessibility analysis of GFP high versus GFP low MAECs identified 501 peaks with reduced accessibility (GFP low peaks) and 3612 peaks with greater accessibility (GFP high peaks) in GFP high MAECs (FDR < 0.05, Figure 2B and Supplementary file 3). For WT and ECKO counterparts, this analysis identified 303 peaks with reduced accessibility (S1pr1 WT peaks) and 472 peaks with enhanced accessibility (S1pr1 ECKO peaks) in S1pr1 ECKO MAECs ( Figure 2B and Supplementary file 3). The~7 fold higher number of GFP high -enriched peaks suggests that GFP high MAECs are more 'activated' (elevated number of chromatin remodeling events) and/or are a heterogeneous mixture of EC subtypes.
To identify relevant transcription factors (TFs), we used the Hypergeometric Optimization of Motif EnRichment (HOMER) (Heinz et al., 2010) suite of tools to reveal over-represented motifs in each set of differentially accessible peaks (DAPs). GFP high peaks were enriched with p65-NFkB, AP-1, STAT3, SOX17, COUP-TFII, and NUR77 motifs ( Figure 2C and Supplementary file 3). In contrast, GFP low peaks showed reduced occurrence of these motifs ( Figure 2D). S1pr1 ECKO peaks were enriched with p65-NFkB motifs, while S1pr1 WT peaks were markedly enriched with glucocorticoid response elements (GREs) and modestly enriched with STAT3, GATA2, ATF1, SOX17 and COUP-TFII motifs ( Figure 2E and F). Examination of ATAC-seq reads centered on selected binding sites (p65, NUR77, COUP-TFII, ATF1, GATA2, and GRE) showed local decreases in accessibility at motif centers, suggestive of chromatin occupancy by these factors ( Figure 2G and H).
We used the ATAC-seq footprinting software HINT-ATAC (Li et al., 2019) to assess genomewide putative chromatin occupancy by TFs. HINT-ATAC identified enhanced footprinting scores at NFKB1 and NFKB2 motifs in S1pr1 ECKO MAECs, whereas GFP high MAECs showed increased scores at RELA motifs and to a lesser extent at NFKB1 and NFKB2 motifs (Figure 2-figure supplement 3A and B). This analysis also identified motifs of the TCF/LEF family (LEF1, TCF7, TCF7L2) as GFP high -enriched, but not S1pr1 ECKO-enriched (Figure 2-figure supplement 3A-D). Consistent with HOMER analysis of DAPs, HINT-ATAC identified COUP-TFII, NUR77, TCF4, and SOX17 motifs as exhibiting greater footprinting scores in GFP high MAECs, while GFP low MAECs showed enhanced putative chromatin occupancy at ATF1 motifs.
Analysis of DAPs showed that only the p65-NFkB motif was commonly enriched between GFP high and S1pr1 ECKO MAECs. This observation is consistent with our RNA-seq analysis, which identified cytokine/NFkB pathway suppression by S1PR1 signaling in MAECs ( Figure 1D). Enrichment of COUP-TFII, NUR77, and AP-1/bZIP motifs in open chromatin of GFP high MAECs, but not S1pr1 ECKO MAECs, further suggests that high levels of S1PR1/ß-arrestin coupling occurs in heterogenous populations of aortic ECs.

Single-cell RNA-seq analysis of GFP high and GFP low MAECs reveals eight distinct EC clusters
Imaging studies demonstrated that GFP high MAECs are restricted to specific anatomical locations. To test the hypothesis that these represent specific EC subpopulations, we employed single-cell RNA-seq (scRNA-seq) on FACS-sorted GFP high and GFP low MAECs. In total, 1152 cells were sequenced (768 GFP high and 384 GFP low ) using the Smart-seq2 protocol. An average of 300,000 aligned reads/cell were obtained and corresponded to~3200 transcripts/cell. Cdh5 transcripts were broadly detected, consistent with endothelial enrichment of sorted cells (Figure 3-figure supplement 1A). S1pr1 and Arrb2 were also broadly detected (Figure 3-figure supplement 1A), suggesting that receptor activation rather than expression of these factors accounts for heterogenous binding sites. ATAC-seq reads were centered on Tn5 cut sites, trimmed to 10 bp, and nucleotide-resolution bigwig files were generated using DeepTools with reads per genomic content (RPGC) normalization. Reads were subsequently centered on TF binding motifs identified in (C-F) and viewed as mean read densities across 600 bp windows. The online version of this article includes the following figure supplement(s) for figure 2:   reporter expression in MAECs. eGFP mRNA was primarily restricted to GFP high MAECs, particularly in the cluster designated aEC1 (Figure 3-figure supplement 1A).
Analysis of GFP high and GFP low MAECs using the velocyto/pagoda2 pipeline (R code in Source Code File 3) (Fan et al., 2016;La Manno et al., 2018) revealed nine clusters upon T-distributed stochastic neighborhood embedding (t-SNE) projection ( Figure 3A). 6 of the nine clusters grouped together in a 'cloud', whereas three clusters formed distinct populations. We used hierarchical differential expression analysis to identify signature marker genes of each cluster ( Figure 3B).
Genes uniquely detected in one of the distinct clusters included vascular smooth muscle cell (VSMC)-specific transcripts such as Myh11, Myom1, and Myocd ( Figure 3B and C; Figure 3-figure supplement 1B). Therefore, this cluster was designated VSMC-like as these cells may represent MAECs sorted along with fragments of VSMCs, or 'doublets' of ECs and VSMCs. Because these cells may represent contamination in an otherwise pure pool of aortic ECs, we omitted this VSMC-like cluster from subsequent analyses. The remaining eight EC clusters were further analyzed.
Lymphatic EC (LEC) markers such as Flt4 (VEGFR3), Prox1, and Lyve1, as well as the venous marker Nr2f2 (COUP-TFII), were detected in a distinct cluster ( Figure 3B  We individually compared LECs, vECs, and VSMCs to the remainder of ECs as a 'pseudo-bulk' cluster to generate a list of transcripts enriched (Z-score >3) for each of these three clusters. We performed the same analysis for aEC1, aEC2, aEC3, aEC4, aEC5, and aEC6, but used only arterial ECs as the comparator. For example, aEC1-enriched transcripts were identified by generating a 'pseudobulk' merge of all aEC2, aEC3, aEC4, aEC5, and aEC6 cells, while aEC2-enriched transcrips were compared to the pseudo-bulk merge of aEC1, aEC3, aEC4, aEC5, and aEC6. The top 32 transcripts that resulted from this analysis are shown in Representative marker genes of aEC1, aEC2, aEC3, and aEC4 are shown in t-SNE embedding in Figure 3G-J.
We note that the abovementioned transcripts ( Figure 3-figure supplement 4B) were among those most differentially expressed between GFP high and GFP low MAECs by bulk RNA-seq (Figure 1-figure supplement 1C). This is demonstrative of consistency between our bulk and singlecell datasets. For functional insights into arterial EC clusters, we analyzed aEC1-aEC4 enriched transcripts with the Gene Set Enrichment Analysis (GSEA) tool ( Figure 4A-D and Supplementary file 5). aEC1 cells were enriched with transcripts associated with GPCR/MAPK signaling (Rasgrp3, Rapgef4, Rgs10, Mapk4k3, S1pr1) as well as VEGF, integrin, and tight-junction pathways (Flt1, Vegfc, Pgf, Igf2, Vcan, Sema3g, S100a4, Jam2, Cldn5). The aEC2 cluster presented a different profile with enriched terms related to immune/inflammatory pathways, TGFß signaling and mRNA processing. Elevated expression of Vcam1, Icam1, Traf6, Cxcl12 and NFkb1 may suggest that these ECs represent an inflammatory cluster.
In contrast, aEC3 cells were enriched with 'immediate-early' transcripts, including those of the AP-1 transcription factor family (Atf3, Jun, Jund, Junb, Fos, Fosb). . Notably, a recent study of young (8 week) and aged (18 month) normal mouse aortic endothelium also identified a cluster of Atf3-positive cells only in young endothelium, as determined by both scRNA-seq and immunostaining for ATF3 (McDonald et al., 2018). Markers of these cells were also identified as top aEC3 markers (e.g. Fosb, Jun, Jund, Junb, Dusp1) suggesting that aEC3 cells are similar to the regenerative Atf3-positive cluster described by McDonald et al. (2018). aEC4 cells were enriched with transcripts related to cell-ECM interations, glycosaminoglycan metabolism, and collagen formation (Pcolce2, Frzb, Spon1, Col4a4, Mfap5, Hyal2). The other two arterial clusters (aEC5 and aEC6) appeared less distinctive (i.e. harbored relatively few enriched transcripts with high Z-scores and fold change values) and therefore were not analyzed. Collectively, these data suggest that scRNA-seq identified more than four distinct MAEC subtypes with unique transcriptomic signatures.
Next, we compared our dataset with information from two recent scRNA-seq studies that also sub-categorized ECs of the normal mouse aorta (Kalluri et (Lukowski et al., 2019) and 'EC1' (Kalluri et al., 2019), strongly resembled the aggregate of aEC2-6 as they were enriched with Gata6, Vcam1, Dcn, Mfap5, Sfrp1, Eln, and Cytl1 (Figure 4-figure supplement 1A and B). Taken together, these information from three independent groups strongly suggest that a major source of heterogeneity in the normal adult mouse aorta, as revealed by scRNA-seq analysis, includes differences between LEC, vEC, aEC1 and aEC2-6.
In contrast, TSP1 staining was exluded from ITGA6-expressing GFP high MAECs at branch points ( Figure 5F; see also Similar to TSP1, VCAM1 staining was low at branch points, while surrounding cells showed heterogeneous levels of immunoreactivity ( Figure 5H). We identified isolated GFP+ cells and patches of GFP+ cells distal from branch points with high VCAM1 immunoreactivity (Figure 5-figure supplement 2A). As expected, intraperitoneal LPS administration induced VCAM1 expression uniformly in aortic ECs ( Figure 5-figure supplement 2B). These data suggest that the aEC1 population includes MAECs that are anatomically specific to branch points, while the aEC2 population is located distally and heterogeneously from branch points but nonetheless exhibits S1PR1/ß-arrestin signaling.

Analysis of branch point-specific arterial ECs cluster aEC1
Having identified the aEC1 cluster as including cells which demarcate thoracic branch points orifices, we sought to characterize these unique cells in further detail. As shown in Figure 4A, aEC1-enriched transcripts are associated with a diverse range of signaling pathways, including MAPK/GPCR, VEGF, and integrin signaling. Among aEC1-enriched transcripts, 16 were up-regulated in S1pr1 ECKO MAECs, five were down-regulated (including S1pr1) and the remaining 390 were not differentially expressed ( Figure 6A; see also Supplementary file 4). The 16 ECKO up-regulated transcripts included positive regulators of angiogenesis (Pgf, Apold1, Itga6, Kdr) (Mirza et al., 2013;Olsson et al., 2006;Primo et al., 2010), regulators of GPCR signaling (Rgs2, Rasgrp3), and Cx3cl1 (Fractalkine), which encodes a potent monocyte chemoattractant (White and Greaves, 2012). We noted that several aEC1-enriched transcripts were also expressed in LEC and vEC, but nonetheless were depleted from the remainder of arterial ECs (aEC2-6) ( Figure 6B). We examined the 25 transcripts most specific to aEC1 (log 2 [fold-change] vs. all ECs > 4) and observed that 5 of these (Dusp26, Eps8l2, Hapln1, Lrmp, and Rasd1) showed differential expression upon loss of S1PR1 function in MAECs ( Figure 6C). Therefore, the majority of transcripts specific to aEC1 do not appear to require S1PR1 signaling for normal levels of expression. For example, protein levels of the aEC1 marker ITGA6 were not markedly affected at branch point orifices in S1pr1 ECKO animals ( Figure 6D). The~2 fold increase in Itga6 transcript levels in S1pr1 ECKO MAECs may be due to expression in heterogeneous (non aEC-1) populations, reflecting increases in LECs and/or aEC2-6. Nonetheless, these data suggest that S1PR1 signaling is required for normal expression of some, but not the majority, of transcripts enriched in aEC1 cells located at orifices of intercostal branch points.
We addressed whether circulatory S1P is required for S1PR1/ß-arrestin coupling in arterial aortic endothelium by genetically deleting the two murine sphingosine kinase enzymes, Sphk1 and Sphk2. We generated S1PR1-GS mice harboring Sphk1 f/f and Sphk2 -/alleles, bred this strain with the tamoxifen-inducible Rosa26-Cre-ER T2 allele, and induced Sphk1 deletion by tamoxifen injection into adult mice (see Materials and methods). 5-6 weeks after tamoxifen administration, plasma S1P  . Gene expression in aEC1 cells is largely independent of S1P/S1PR1 signaling. (A) Pie chart of all aEC1enriched transcripts (versus aEC2-6) indicating those which were also differentially expressed in S1pr1 ECKO MAECs. (B) Heatmap (row Z-scores) of the 20 transcripts that were both S1PR1-regulated and aEC1-enriched (left). Expression of these transcripts in each cluster from scRNA-seq analysis is shown (right). (C) Gene expression in aEC1 was compared to all ECs (LEC, vEC, and aEC2-6 collectively). All transcripts expressed greater than 16-fold higher in aEC1 are shown (25 transcripts total). Red, blue and black transcript names indicate up-regulated, downregulated or similar levels of expression in S1pr1 ECKO MAECs, respectively. (D) Immunostaining of S1pr1 ECKO and WT aortae whole-mount en face preparations for ITGA6 and VEC. Images are representative of observations from 2 pairs of animals (N = 2). (E) ITGA6 and VEC immunostaining of whole-mount en face preparations of Figure 6 continued on next page concentrations in Cre-animals were 631 ± 280 nM, whereas S1P was undetectable in plasma from Cre+ mice. Cre+ mice had~7 fold fewer non-branch point GFP+ ECs ( Figure 6E and F), suggesting that S1PR1/ß-arrestin coupling in aEC2 cells is ligand-dependent. In contrast, the number of branch point GFP+ EC was not significantly different between Cre+ and Cre-animals, suggesting ligandindependent S1PR1/ß-arrestin coupling in aEC1 population ( Figure 6E and F). Furthermore, ITGA6 expression at branch point orifices was unaltered in Cre+ mice and therefore is independent of circulatory S1P. Taken together, these data suggest that the unique transcriptome of aEC1 cells is largely independent of S1P/S1PR1 signaling while in aEC2 (and perhaps others), ligand-dependent S1PR1 signaling predominates.
Spatio-temporal and molecular differences between branch pointspecific arterial EC cluster aEC1 and non-branch point ECs . Transcripts encoding other TFs, such as Tcf4, Ets1, Sox18, Epas1, Mef2c, and Tox2, were aEC1-enriched (Z-score >3 and<7) but more heterogeneously distributed in other arterial clusters ( Figure 7A and Figure 7-figure supplement 1A). These data are consistent with our ATAC-seq analysis, which showed over-representation of SOX17, TCF4, and NUR77 motifs in chromatin specifically open in GFP high MAECs. In contrast, Gata6 (ubiquitous among aEC2-6) and Gata3 (heterogeneous among aEC2-6) were both depleted from aEC1 ( Figure 7A and Immunostaining of thoracic aorta en face preparations for LEF1 showed nuclear immunoreactivity in GFP high ECs at branch point orifices but not in adjacent ECs ( Figure 7B). We noted that all LEF1+ cells also exhibited ITGA6 expression, confirming these two proteins as markers of aEC1 cells at branch point orifices ( Figure 7B). These data suggest involvement of LEF1, a downstream TF of Wnt/ß-catenin signaling, in regulating gene expression in aEC1 cells.
For broader insight into TF activity near aEC1 genes, we extracted all GFP high and GFP low merged peaks that intersected a 100 kb window centered on the TSSs of all aEC1-enriched (Z-score >3) and -depleted (Z-score < À3) genes. HOMER analysis of these peaks showed enrichment of SOX17, 'Ets1-distal', MEF2C, and NUR77 motifs near the aEC1-enriched genes (Figure 7-figure supplement 1C). In contrast, GATA motifs (GATA2, GATA3, GATA6), as well as the NKX2.2 motif, were enriched near aEC1-depleted genes (Figure 7-figure supplement 1C). Collectively, these data suggest roles for specific transcription factors, such as LEF1, SOX17, NUR77, and GATA6, in mediating transcriptional events which distinguish aEC1 from the remainder of aortic arterial ECs.
Examination of postnatal day 6 (P6) S1PR1-GS aortae showed that ECs at branch point orifices expressed GFP and ITGA6 in a manner similar to adult S1PR1-GS mice ( Figure 7C). We noted that non-branch point GFP+ EC were less frequent in P6 mice relative to young adults ( Figure 7C). Quantification of non-branch point GFP+ EC in aortae of P6, P12, young adult (3-4 months), and 15 month old mice revealed a~5 fold increase from P12 to young adult ( Figure 7D and There was no appreciable difference in non-branch point GFP+ EC frequency between P6 and P12 or between young adult and 15 months. In contrast, the number of GFP+ EC at branch points was stable from P12 throughout adulthood (Figure 7-figure supplement 2E). These data suggest that S1PR1/ß-arrestin coupling in non-branch point ECs increases with post-natal age. thoracic aortae from S1PR1-GS mice bearing Sphk2 -/-Sphk1 f/f or Sphk2 -/-Sphk1 f/f Rosa26-Cre-ER T2 alleles. (F) Quantification and statistical analyses (unpaired t-test) of GFP+ EC at branch point (12 branches) and non-branch point (six fields) locations from 2 pairs of mice (N = 2). Quantified fields are as shown in (E). Branch point EC were defined as cells included in the first three rows around the edge of orifices. Only GFP+ EC in the same Z-plane as surrounding arterial EC were counted (GFP+ EC of intercostal arteries were in a different Z-plane and therefore were not counted as branch point EC). Scale bars are 100 mM. To test the hypothesis that aEC2-like (GFP high ) cells occur at greater frequency during aging, aortae from P6, young adult, and 15 month old mice were immunostained for VEC and TSP1, which is abundantly expressed by aEC2 cells. TSP1 immunoreactivity was markedly reduced to near-undetectable levels in P6 mice relative to older counterparts ( Figure 7E and F). Taken together, these data suggest that TSP1-expressing (aEC2-like) cells increase in frequency while aEC3-like cells disappear over time in the aorta intima, which was described previously (McDonald et al., 2018).
We sought to determine if aEC1-enriched transcripts other than Itga6 and Thbs1 show evidence of temporally-regulated expression in aortic endothelium. Examination of the embryonic gene expression database (Cao et al., 2019) revealed that other aEC1-enriched transcripts, such as Igf2, Flt1, Slc6a6, and Lef1 were expressed at greater levels in embryonic ECs relative to aEC1-depleted transcripts (e.g. Vcam1, Frzb, Pcolce2, Cyp1b1, Sod3) (Figure 7-figure supplement 3D) that are enriched in aEC2 and/or aEC4. Thus, aortic endothelium exhibits spatio-temporal regulation of genes such as Itga6 and Thbs1, with expression of aEC1 genes such as Itga6 becoming restricted to branch point orifices over time while expression of aEC2-enriched genes (e.g. Thbs1) increases throughout the aorta intima after P6.
We further explored potential roles for LEF1, SOX17, NUR77, and GATA6 in aEC1 gene expression by examining binding sites for these factors near genes encoding aEC1-enriched transcripts. A prominent GFP high -specific peak in the first intron of Itga6 harbored a LEF1 motif ( Figure 7I). In followed by Bonferroni's post hoc test, * P = 0.04; ** P = 0.04; *** P = 0.0008. (E) Representative images of TSP1 and VEC immunostaining of thoracic aorta en face preparations from P6, 3 month old, 15 month old mice. Scale bars are 200 mM. (F) Quantification of TSP1 pixels normalized to total VEC pixels per field from P6 (n = 3), 3 month old (n = 2), and 15 month old (n = 2) mice. One-way ANOVA followed by Bonferroni's post hoc test, ** P = 0.001; **** P < 0.0001. (G) Sagittal cryosection (14 mM) of an E11.5 S1PR1-GS mouse embryo immunostained for CD31, ITGA6, and TSP1. Arrows indicate the dorsal aorta. Scale bar is 500 mM. (H) The dorsal aorta at higher magnification with CD31, ITGA6, and GFP channels shown. Scale bar is 100 mM. (G) and (H) are representative of observations from two embryos. (I-J) Genome browser image of GFP high and GFP low MAECs ATAC-seq signal (RPGC normalized) at Itga6 and Thbs1 (TSP1) loci. Peaks with increased accessibility in GFP high MAECs (GFP high DAPs) are indicated in (I). All GFP high and GFP low MAECs peaks are shown in (J). All LEF1, nuclear receptor subfamily four group A member 1 (NUR77), transcription factor SOX-17 (SOX17), and transcription factor GATA-6 (GATA6) motifs identified in consensus peaks (Figure 2A) are also shown. Orange bars in (I) highlight GFP high MAECs DAPs containing LEF1, SOX17, and/or NUR77 motifs. Yellow bars in (J) highlight peaks containing GATA6 motifs. The online version of this article includes the following figure supplement(s) for figure 7:     addition, a GFP high -enriched peak~115 kb upstream of the transcription start site (TSS) harbored a SOX17 motif. The Nog gene contained SOX17 and NUR77 motifs in upstream GFP high -specific peaks (Figure 7-figure supplement 4A). Each of these putative enhancers lacked GATA6 motifs. In contrast, open chromatin near the Thbs1 (TSP1) gene lacked SOX17 and LEF1 sites and instead contained three GATA6 sites ( Figure 7J). Furthermore, the Frzb and Pcolce2 genes, which encode aEC4-enriched transcripts and were up-regulated in GFP low MAECs (Figure 1-figure supplement  2D), harbored intronic GFP low MAEC-specific peaks containing GATA6 binding sites (Figure 7-figure supplement 4B).
Together, these findings characterize the aortic branch point-specific arterial EC subpopulation designated aEC1. The data strongly suggest that the GFP high status of MAECs at branch point orifices, as well as their unique gene expression program, is not dynamic throughout postnatal life. Rather, the gene expression specification of these cells likely occurs during development in tandem with epigenetic changes (i.e. chromatin accessibility) and is stable throughout adulthood. These cells have a unique anatomical location in postnatal mice and exhibit high S1PR1/ß-arrestin coupling. The transcriptome of this EC subpopulation does not appear to be directly regulated by S1PR1 signaling. Rather, a combination of TFs, such as NUR77, NURR1, SOX17, HEY1, and LEF1, likely regulate cluster-defining transcripts in these cells.

S1PR1 signaling in adventital lymphatic ECs regulates immune and inflammatory gene expression
To locate the aorta-associated LEC with a high frequency (97%) of S1PR1/ß-arrestin coupling, we utilized antibodies against LYVE1 and VEGFR3 (Flt4). LYVE1 marks most but not all LEC subtypes and VEGFR3 is a pan-LEC marker (Wang et al., 2017). Sagittal sections of the S1PR1-GS mouse thoracic aorta revealed that a subset of adventitial LYVE1+ LECs are GFP+ and thus exhibit S1PR1/ß-arrestin coupling ( Figure 8A).
In addition to expression in aEC1, Itga6 transcripts were detected in vEC and LEC populations (Figure 8-figure supplement 1A). Consistently, we observed GFP+ITGA6+LYVE1+ LEC on the adventitial side of the aorta, in proximity to a GFP+ITGA6+LYVE1-arterial branch point (Figure 8figure supplement 1C).
We examined the role of LEC-derived S1P in LEC S1PR1/ß-arrestin coupling by generating S1PR1-GS mice deficient in lymph S1P (S1pr1 ki/+ :Sphk1 f/f :Sphk2 -/-:Lyve1-Cre + ) (Pham et al., 2010). LECs were identified in mesenteric vessels by immunostaining for the LEC-specific transcription factor prospero homeobox protein 1 (PROX1). Collecting lymphatics, which exhibit coverage of asmooth muscle actin (ASMA) positive cells (Wang et al., 2017), were identified by ASMA co-staining. H2B-GFP control mice showed low levels of GFP expression in lymphatic valves and the remainder of LECs were GFP-( Figure 8C). In stark contrast, GFP+ LECs were observed at high frequency (73%) in ASMA+ structures but not in ASMA-structures (6%) of S1PR1-GS mice ( Figure 8D and E). This suggests that S1PR1/ß-arrestin coupling occurs primarily in collecting lymphatics. Mice lacking S1P in lymph (Lyve1-Cre + ) exhibited a 10-fold reduction in GFP+ LEC in ASMA+ structures (7%), indicating that the ligand S1P induces S1PR1/ß-arrestin coupling in collecting lymphatic vessels. These data are consistent with abundant LEC expression of the S1P transporter Spns2 (Figure 3figure supplement 4A), which is also required for normal levels of lymph S1P (Simmons et al., 2019). Taken together, scRNA-seq analysis of aortic ECs identified two anatomically distinct arterial EC populations (branch point and non-branch point), as well as a collecting lymphatic-like adventitial LEC population, each of which shows high S1PR1/ß-arrestin coupling.
For insight into S1PR1-mediated gene expression in aorta-associated LECs, we divided S1pr1 ECKO up-regulated genes ( Figure 1C and Figure 1-figure supplement 2E) according to their cluster assignments from scRNA-seq analysis. 30% of up-regulated genes were enriched in the LEC Figure 8. S1PR1/ß-arrestin coupling in aorta-associated lymphatic vessels and S1P-dependent signaling in LECs. (A) Representative confocal image of a sagittal cryosection (14 mM) of a S1RP1-GS mouse aorta immunostained for VEC and lymphatic vessel endothelial hyaluronic acid receptor 1 (LYVE1) (N = 2). White arrows indicate GFP+ arterial ECs at a branch point orifices and yellow arrows indicate adventitia-associated GFP+ LECs. Scale bar is 50 mM. (B) Confocal image of a S1RP1-GS mouse thoracic aorta whole-mount preparation immunostained for vascular endothelial growth factor receptor 3 (VEGFR3; Flt4), LYVE1, and CLDN5 with the tunica intima in contact with coverslip. Orange stars indicate VEGFR3+LYVE1+GFP+ areas, cyan stars indicate VEGFR3+LYVE1 low GFP+ areas, Figure 8 continued on next page cluster (Z-score >3 versus the remainder of ECs), while only 7% were enriched in vEC and/or aEC1-6 ( Figure 9A and Suplementary File 6). None of the S1pr1 ECKO down-regulated transcripts were enriched in the LEC cluster (see Supplementary file 6). A heatmap of the 78 LEC transcripts up-regulated in S1pr1 ECKO MAECs is shown in Figure 9B. Among these were chemokine/cytokine pathway genes (Irf8, Lbp, Il7, Il33 Ccl21, Tnfaip8l1) as well lymphangiogenesis-associated genes (Kdr, Prox1, Lyve1, Nr2f2) ( Figure 9B and Figure 9-figure supplement 1A). This suggests that loss of S1PR1 signaling in LECs alters transcriptional programs associated with lymphangiogenesis and inflammation/immunity.
Lymphatic vessels associated with the aorta and large arteries, although not well studied, are thought to be involved in key physiological and pathological processes in vascular and immune systems (Csányi and Singla, 2019; Galkina and Ley, 2009). For example, LEC expression of CCL21 mediates dendritic cell recruitment to lymphatic vessels during homeostasis and pathological conditions (Vaahtomeri et al., 2017). This chemokine was abundantly expressed in the LEC cluster (Figure 8-figure supplement 1A) and was up-regulated by S1PR1 loss-of-function ( Figure 9B and Figure 9-figure supplement 1A). Whole-mount staining of S1PR1-GS mouse thoracic aortae for LYVE1 and CCL21 revealed aorta-associated CCL21+ lymphatics with high levels of S1PR1/ß-arrestin coupling ( Figure 9C and D). We noted that CCL21 protein appeared as peri-nuclear puncta that likely marks the trans-Golgi network, as observed in dermal LECs (Vaahtomeri et al., 2017). Sagittal sections of the thoracic aorta indicate that GFP high , LYVE1 + adventitial lymphatics express CCL21 protein ( Figure 9D), suggesting that ß-arrestin-dependent down-regulation of S1PR1 correlates with CCL21 expression ( Figure 9D).
These studies show that S1PR1/ß-arrestin coupling in a subset of adventitial lymphatics correlates with S1PR1-mediated attenuation of lymphagiogenic/inflammatory gene expression. We observed that the fraction of PDPN + LECs was not altered between S1pr1 ECKO and WT aorta tissues (Figure 9-figure supplement 1B), indicating that there is not widespread lymphagiogenesis or LEC proliferation. However, analysis of Flt4 vs. Pdpn expression and Ccl21a vs. Pdpn expression revealed the presence of LECs expressing Flt4 and Ccl21a but not Pdpn (Figure 9-figure supplement 2). Therefore, it is possible that S1pr1 ECKO leads to expansion of a PDPN low population in the adventitia. Nonetheless, these data indicate that S1PR1 mediates gene expression in adventitial lymphatics of the adult mouse aorta under homeostasis.

Discussion
Intracellular signaling through G protein-and ß-arrestin-dependent pathways is tightly regulated at the levels of GPCR expression and ligand availability. Endothelium of major organs, such as brain, lung, skeletal muscle, and the aorta, express unique sets of GPCRs with only five receptors commonly expressed, one of which is S1pr1 (Kaur et al., 2017). Despite ubiquitoius endothelial S1pr1 expression, S1PR1/ß-arrestin coupling in vivo, as reported by GFP in S1PR1-GS mice, revealed heterogeneous signaling in multiple organs (Kono et al., 2014;Galvani et al., 2015). Here, we describe high frequency of aortic GFP+ (i.e. S1PR1/ß-arrestin coupling) ECs around intercostal branch point orifices and heterogeneous GFP+ ECs throughout the remainder of intimal aortic endothelium. We and others have shown that endothelial ablation of S1pr1 (S1pr1 ECKO) or reduction of circulatory and magenta stars indicate VEGFR3+LYVE1+GFP-areas. Image is representative of observations from three mice (see also Figure 8-figure supplement 1B). Scale bar is 100 mM. (C-D) Representative confocal images of mesenteric lymphatics from H2B-GFP control (C) and S1PR1-GS mice bearing Sphk2 -/-:Sphk1 f/f (Cre-) or Sphk2 -/-: Sphk1 f/f :Lyve1-Cre (Lyve1-Cre+) alleles whole-mounted and immunostained for PROX1 and ASMA. White arrows indicate ASMA+PROX1+GFP+, yellow arrows indicate PROX1+GFP-, cyan arrows indicate ASMA+PROX1+GFP-. Scale bars are 100 mM. (E) Quantification of GFP+ signal over ASMA+PROX1+ and ASMA-PROX1+ areas from Cre-(n = 5), Lyve1-Cre (n = 4), and H2B-GFP control (n = 2) mice. One-way ANOVA followed by Bonferroni's post hoc test, **** P < 0.0001. The online version of this article includes the following figure supplement(s) for figure 8:  Confocal images of two fields (c' and c'') of a whole-mounted S1RP1-GS mouse thoracic aorta immunostained for LYVE1 and C-C motif chemokine 21 (CCL21). The tunica adventitia was facing the coverslip. (D) Confocal images of a sagittal cryosection (14 mM) of an S1RP1-GS mouse thoracic aorta immunostained CCL21 and LYVE1. Yellow Figure 9 continued on next page S1P disrupts endothelial barriers (Camerer et al., 2009;Christensen et al., 2016;Christoffersen et al., 2011;Oo et al., 2011;Yanagida et al., 2017). Furthermore, the descending aorta of S1pr1 ECKO mice displayed exacerbated plaque formation in the Apoe -/-Western diet (WD)-induced atherosclerosis model (Galvani et al., 2015). However, we lack information regarding EC transcriptional responses that correlate with or are directly downstream of S1PR1 signaling. Here, we profiled the transcriptomes and open chromatin landscapes of aortic ECs with high (GFP high ) or low (GFP low ) levels of S1PR1/ß-arrestin coupling, as well as S1pr1 ECKO aortic ECs. S1pr1 ECKO MAECs up-regulated transcripts in TNFa/cytokine signaling pathways and exhibited enhanced chromatin accessibility at NFkB binding sites. Concomitantly, the glucocorticoid receptor pathway was suppressed. These mRNA and chromatin signatures were shared between S1pr1 ECKO and GFP high MAECs, suggesting that persistent ß-arrestin recruitment to S1PR1 can result in downregulation of membrane-localized S1PR1 and a subsequent loss-of-function phenotype.
There were many (2,145) DEGs between GFP high and GFP low MAECs, but relatively few (365) between ECKO and WT MAECs, which suggested that the GFP high and/or GFP low populations reported aortic EC subtypes in addition to S1PR1-regulated transcripts. Indeed, chromatin regions uniquely open in GFP high MAECs were enriched with binding sites for TFs that have well-defined but divergent roles in endothelial cells, such as SOX17 (arterial) (Corada et al., 2013;Zhou et al., 2015) and COUP-TFII (venous/lymphatic) (Lindskog et al., 2014). Nonetheless, we present the first collection of putative regulatory regions from freshly isolated mouse aortic ECs, which is a critical dataset for future studies of individual enhancer functionalities.
Our scRNA-seq analysis addressed with high resolution the heterogeneity among MAECs. We identified six arterial EC clusters (aEC1-6), one lymphatic EC cluster (LEC) and one venous EC cluster (vEC). Immunohistochemical analyses revealed LEC cells as including lymphatic structures of the aortic adventita, aEC1 cells as circumscribing intercostal branch point orifices, and aEC2 cells as heterogeneously dispersed throughout intimal endothelium. Each of these clusters exhibited a high frequency (>90%) of GFP high EC. We also describe aEC3 cells, which contained comparatively few (<60%) GFP high ECs. aEC3 cells strongly resembled an Atf3-positive cluster reported by McDonald et al. (2018) that mediates endothelial regeneration (McDonald et al., 2018). Considering that Atf3-positive ECs were absent in old (18 month) mice (McDonald et al., 2018), we hypothesize that aEC3-like cells disappear over time while aEC2-like cells increase in frequency in the aorta intima. This notion is supported by the higher frequency of non-branch point GFP high and TSP1expressing intimal ECs in adult relative to P6 and P12 mice.
When compared with findings from two recent aorta scRNA-seq studies (Kalluri et al., 2019;Lukowski et al., 2019), our clustering segregated vEC, LEC, and aEC1 from aEC2-6. We suspect that use of S1PR1-GS mice facilitated deconvolution of LEC and vEC from the distinct aEC1 population. Despite their proximity to intercostal branch points, aEC1 cells do not exhibit a transcriptomic signature prototypical of inflammation, as might be expected of ECs in an environment with disturbed flow (Chiu and Chien, 2011). Furthermore, high levels of S1PR1/ß-arrestin coupling and expression of unique genes (e.g. Itga6) in aEC1 cells are independent of both circulatory S1P and age in postnatal mice. To further explore temporal regulation of cluster-specific genes, we examined scRNA-seq of FACS-sorted VEGFR2 high cells from E8.25 embryos (Pijuan-Sala et al., 2019) and the endothelial cluster from scRNA-seq of E9.5 to E13.5 embryos (Cao et al., 2019). These embryonic cells exhibited expression of aEC1-enriched transcripts (Lef1, Itga6, Alpl, Flt1, Igf2), but depletion of aEC2-6-enriched transcripts (Thbs1, Sod3, Vcam1, Sfrp1, Pcolce2, Dcn). Together with embryonic EC gene expression data (Figure 7-figure supplement 3) and immunohistochemical analysis of E11.5 S1PR1-GS embryos, our data suggest that aEC1 cells are more characteristic of embryonic EC than are the majority of intimal ECs. Similarly, vEC/LEC-specific transcripts (Prox1, Lyve1, Nr2f2) were also expressed in embryonic ECs in both studies (Pijuan-Sala et al., 2019;Cao et al., 2019). Thus, transcriptomic similarities between aEC1 and LEC/vEC in the adult aorta may be retained from development, perhaps through epigenetic modifications common between these cell types. It is also possible that the unique anatomical location of aEC1 (at the circumferential ridge at aortic branch points) may promote a distinct EC phenotype either because of spatial positioning or environmental factors. S1pr1 ECKO and S1pr1 -/animals display an aortic hyper-branching phenotype between E11.5 and E13.5 that is incompatible with life after E14.5 (Gaengel et al., 2012). Therefore, S1PR1 is required for normal embryonic branching morphogenesis. Consistently, E9.5-E10.5 S1PR1-GS embryos show high GFP expression (S1PR1/ß-arrestin coupling) in the dorsal aorta (Kono et al., 2014). This contrasts with adult aortae, wherein the highest levels of S1PR1/ß-arrestin coupling are concentrated around the orifices of intercostal branch points. These findings further suggest that the unique gene expression program of aEC1 cells is established during morphogenesis of intercostal arteries during development. It is possible that the chromatin accessibility differences that we observed near aEC1-enriched and -depleted genes were established during embryonic or early postnatal life and were retained into adulthood. While we cannot exclude shear forces as contributing to the aEC1 transcriptional program, our data suggest that disturbed flow is not the main driver of gene expression in these cells.
The extent to which aEC1 cells at branch point orifices are functionally distinct remains to be determined. Identification of ITGA6 as a marker of this population will facilitate future studies. For example, combinations of pan-EC, LEC, and ITGA6 antibodies can be used to purify or enrich this population in developmental or disease models (e.g. atherosclerosis). Moreover, cis-elements proximal to aEC1-specific genes can be applied to the 'Dre-rox/Cre-loxP' system (Pu et al., 2018) to specifically manipulate gene expression in aEC1 cells.
Non-branch point (i.e. aEC2) cells require circulatory S1P for S1PR1/ß-arrestin coupling. Similarly, mesenteric LECs required S1P in lymph for normal levels of S1PR1/ß-arrestin coupling. Notably, we did not detect S1pr1 ECKO down-regulated transcripts in LEC cluster, but we did detect five such transcripts in aEC1 cells (Dusp26, Enah, Eps8l2, Hapln1 and S1pr1) and two transcripts in aEC2 cells (Nfix, Znrf2). Among the clusters identified in this study, LEC-enriched transcripts were affected the most upon EC ablation of S1pr1 (Suplementary File 6). This suggests that loss of S1P/S1PR1 signaling either alters cell-intrisic phenotypes of peri-aortic LECs or induces expansion of one or multiple LEC subtypes.
There is accumulating direct and indirect evidence for key roles of adventital lymphatics in atherogenesis (Csányi and Singla, 2019; Maiellaro and Taylor, 2007). For example, auto-antibodies against oxidized LDL (OxLDL) inhibit macrophage OxLDL uptake and mitigate atherosclerosis (Shaw et al., 2000). This implies that antigen presenting cells (APCs) phagocytose OxLDL epitopes, then travel via adventitial lymphatics to lymphoid organs (e.g. lymph nodes) and present OxLDL antigens to B cells. Murine atherosclerotic lesions were found to harbor 'atypical, lymphatic-like' capillaries that were VEGFR3+ but LYVE1- (Taher et al., 2016), which is consistent with our observations of adventitial LEC heterogeneity. Considering the critical role of lymphatic EC-derived CCL21 in regulating the trafficking of APCs (Vaahtomeri et al., 2017), and perhaps other adaptive immune cells, there is an impetus to determine the extent to which adventitial lymphatics are a viable target for atherosclerosis therapy.
Identification of embryonic LEC with high S1PR1/ß-arrestin coupling may suggest a functional role for S1PR1 in developmental lymphangiogenesis. We speculate that future studies of S1pr1 f/f/ mice bearing LEC-specific Cre-drivers will provide insight into S1PR1-mediated events during developmental lymphangiogenesis.
While there is scant information about the roles of S1P/S1PR1 signaling in adult lymphatic vasculature, our findings lay a groundwork for future studies of S1PR1-mediated LEC phenotype regulation in homeostatic processes and inflammatory/autoimmune diseases. A recent study found that lymph-derived S1P facilitates CCL21 deposition in high endothelial venules and dendritic cell recruitment (Simmons et al., 2019). While S1pr1 ECKO animals exhibit exacerbated atherosclerosis (Galvani et al., 2015), we cannot discern whether this was due to phenotypes of lymphatic ECs, arterial ECs, or both cell types. Future studies should use artery-and lymphatic-specific Cre-drivers to distinguish between the roles of S1PR1 in different types of vasculature. Such mechanistic studies will help to determine the utility of S1PR1 modulators in treating lymphatic-mediated vasculopathies.   Continued on next page is described in Takeda et al. (2007). Tamoxifen was administered to S1PR1-GS-S1P-less mice and Cre-littermate controls as described above. Experiments were performed between 23 and 25 weeks after the final tamoxifen dose. Young adult (aged 8 to 12 weeks) males and females were used for sequencing experiments. Males and females of similar age (7 to 18 weeks) were used for imaging studies, unless indicated otherwise. For examination of VCAM1 in Figure 5-figure supplement 2B, 200 mL lipopolysaccharide (Sigma-Aldrich, L2630) in PBS was injected i.p. (5.5 mg/kg) followed by euthanasia and tissue harvest after 9 hr. For timed matings, embryonic day (E) 0.5 was defined as noon on the date of the vaginal plug and embryos were harvested at E11.5.

Lymphocyte counts
Blood was drawn from the retroorbital venous plexus into EDTA tubes and blood cells enumerated with a Hemavet 950 cell counter (Drew Scientific).

Lung vascular leakage assay
Lung vascular integrity was assessed by administration of 6 mL/g 0.5% Evans Blue dye (Sigma # E2129) via the retro-orbital venous plexus. Two hours later, lungs were perfused via the right ventricle with 10 mL heparin-DPBS, harvested, and Evans Blue was extracted overnight at 55˚C in formamide. Lung accumulation of Evans Blue was quantified by measuring absorbance of the lung extract at 620 nm and expressed as corrected for absorbance at 740 nm.

FACS isolation and single cell sequencing of mouse aortic endothelial cells
After CO 2 euthanasia, the right atrium was opened and the left ventricle was perfused with 10 mL phosphate-buffered saline (PBS) (Corning). Aortae were dissected from the root to below the common iliac bifurcation and transferred into ice ice-cold 1x HBSS (Sigma-Aldrich, H1641). Whole aortae were incubated in HBSS containing elastase (4.6 U/mL, LS002292, Worthington), dispase II (1.3 U/ mL, Roche), and hyaluronidase (50.5 U/mL, Sigma-Aldrich, H3506) at 37˚C for 10 min in wells of a 6well plate. Aortae were then transferred to a 100 mm dish with 1 mL HBSS and minced using small scissors. Minced aortae were transferred to a low protein binding 5 mL tube containing Liberase (0.6 U/mL, Sigma-Aldrich), collagenase II (86.7 U/mL, LS004174, Worthington), and DNase (62.0 U/mL, Sigma-Aldrich, D4527) in 4.3 mL HBSS and incubated at 37˚C for 40 min with rotation in a hybridization oven. The cell suspension was then triturated 10 times through an 18 G needle to dissociate clumps, followed by addition of 400 mL STOP solution (3 mM EDTA, 0.5% fatty acid-free bovin serum albumin (FAF-BSA) (Sigma-Aldrich, A6003) in 1x HBSS). For the remainder of the procedure, cells were kept on ice and all centrifugation steps were performed at 4˚C.
Cells were spun at 500 xg for 5 min, the supernatant was removed, then cells were washed with 4 mL STOP solution and spun at 500 xg for 5 min. The supernatant was removed, then cells were washed with 4 mL blocking solution (0.25% FAF-BSA in HBSS) and filtered through FACS tubes with filter caps (Falcon). After centrifugation and supernatant removal, cells were stained with phycoerythrin (PE)-conjugated anti-mouse CD31 (MEC13.3, Biolegend, San Diego, CA), allophycocyanin (APC)-conjugated anti-mouse CD45 (30-F11, Biolegend) and APC-conjugated TER119 (Biolegend, 116212) antibodies in blocking solution with anti-CD16/32 (2.5 mg/mL) for 45 min on ice. DAPI (0.7 mM) was added for the final 5 min of staining to exclude dead cells. Aortic cells were washed with 4.5 mL FACS buffer (0.25% FAF-BSA in PBS) before sorting for selection of CD31 + /CD45 -/TER119 -/ GFP high and CD31 + /CD45 -/TER119 -/GFP low cells using BD FACSAria II (BD Bioscience) (see Figure 1B). Cells from S1pr1 WT and -ECKO mice were sorted using the GFP low gate (see Figure 1figure supplement 1A) because it includes MAECs from mice genetically negative for the H2B-GFP reporter allele and stained with the same antibody panel. Cells from 2 to 4 aortae of age and sexmatched adult mice were pooled for each individual experiment (ATAC-seq, RNA-seq and scRNAseq). Cells were sorted into either 0.1% FAF-BSA/PBS or buffer RLT (Qiagen) supplemented with bmercaptoethanol for ATAC-seq and RNA-seq, respectively.

Bulk RNA-seq and analysis
Cells sorted into buffer RLT were subjected to total RNA extraction using the RNeasy Micro Kit (Qiagen). The High Sensitivity RNA ScreenTape (Agilent) was used to verify RNA quality before synthesis of double-stranded cDNA from 5 to 10 ng RNA using the SMART-Seq2 v4 Ultra Low RNA Kit for Sequencing (Takara Bio) according to the manufacturer's instructions. Agilent 2100 Bioanalyzer and High Sensitivity DNA Kit (Agilent) were used to verify cDNA quality. cDNA libraries were prepared for sequencing using the Illumina Nextera XT2 kit (Illumina), and~20-40 million paired-end reads (2 Â 75 bp) were sequenced for each sample.

ATAC-seq and analysis
ATAC-seq libraries were prepared according to the previously described fast-ATAC protocol (Corces et al., 2016). Briefly, 800-4,000 FACS-isolated cells in 0.1% FAF-BSA/PBS were pelleted by centrifugation at 400 Âg at 4˚C for 5 min. Supernatant was carefully removed to leave the cell pellet undisturbed, then cells were washed once with 1 mL ice-cold PBS. The transposition mix [25 mL buffer TD, 2.5 mL TDE1 (both from Illumina FC-121-1030), 1 mL of 0.5% digitonin (Promega, G9441) and 16 ml nuclease-free water] was prepared and mixed by pipetting, then added to the cell pellet. Pellets were disrupted by gently flicking the tubes, followed by incubation at 37˚C for 30 min in an Eppendorf ThermoMixer with constant agitation at 300 rpm. Tagmented DNA was purified using the MinElute Reaction Cleanup Kit (Qiagen, 28204), and subjected to cycle-limiting PCR as previously described (Buenrostro et al., 2013). Transposed fragments were purified using the MinElute PCR Purification Kit (Qiagen, 28004) and Agilent DNA Tapestation D1000 High Sensitivity chips (Agilent) were used to quantify libraries.~20-60 million paired-end reads (2 Â 75 bp) were sequenced for each sample on a NextSeq instrument (Illumina).
Paired-end reads were separated, centered on Tn5 cut sites, and trimmed to 10 bp using a custom in-house script (Source Code File 1). Peaks were called using the MACS2 'callpeak' script  with options: -B -keep-dup all -nomodel -nolambda -shift À75 -extsize 150. Reads mapping to murine blacklisted regions and mitochondrial DNA were masked out of peak lists using the Bedtools 'intersect -v' script (Quinlan and Hall, 2010).
Replicates from each biological group were merged using the bedops 'merge' script to generate one high-confidence peak set for each of the four biological groups (GFP high , GFP low , S1pr1 ECKO, S1pr1 WT) (Neph et al., 2012). These four peak sets were then merged to generate a merged, consensus peak set of 123,473 peaks. For each replicate, reads covering consensus peak intervals were counted using the Bedtools 'coverage' script with the '-counts' option (Quinlan and Hall, 2010). The resultant count table was input to edgeR (Robinson et al., 2010) to determine differentially accessible peaks (DAPs).
DAPs were used as input for the HOMER 'findMotifsGenome.pl' script with the option '-size given' to identify motifs enriched in peaks with enhanced accessibility in either GFP high , GFP low , S1pr1 ECKO, or S1pr1 WT MAECs (Heinz et al., 2010).
Nucleotide-resolution coverage (bigWig) tracks were generated by first combining trimmed reads from each replicate, then inputting the resultant. bam files to the DeepTools (Ramírez et al., 2016) 'bamCoverage' script with options '-normalizeUsing RPGC -binSize 1'. Heatmaps of ATAC-seq reads within 600 bp of p65, NUR77, COUP-TFII, ATF1, GATA2, and GRE motifs were generated by centering coverage tracks on each motif identified in DAPs. These motifs were identified using the HOMER script 'annotatePeaks.pl' with the '-m -mbed' options. All heatmaps were generated using DeepTools and all genome browser images were captured using Integrative Genomics Viewer (Robinson et al., 2011). scRNA-seq analysis 1152 Fastq files (one per cell) were aligned to the GRCm38 (mm10) genome assembly using STAR with options -runThreadN 4 -outSAMstrandField intronMotif -twopassmode Basic. Bam files were input to the velocyto (http://velocyto.org/velocyto.py/) command-line script 'run-smartseq2' (La Manno et al., 2018). Expressed repetitive elments were downloaded from the UCSC genome browser and masked from analysis using the '-m' option of the 'run-smartseq2' script. The resultant table of read counts per transcript ('loom' file) was input to the PAGODA2 (https://github.com/hmsdbmi/pagoda2) R package for further analysis (Fan et al., 2016;La Manno et al., 2018). The details of our R code are provided in Source Code File 3.
After variance normalization, the top 3000 overdispersed genes were used for principal component analysis (PCA). An approximate k-nearest neighbor graph (k = 30) based on a cosine distance of the top 100 principal components was used for clustering. Clusters were determined using the multilevel community detection algorithm. PCA results were plotted using the 'tSNE' embedding option of the PAGODA2 'r$getEmbedding' function. Heatmaps of gene expression embedded on hierarchical clustering, differential expression analyses, and expression of individual transcripts on the tSNE embedding were generated using the graphical user interface at http://pklab.med.harvard. edu/nikolas/pagoda2/frontend/current/pagodaLocal/. We generated the binary (.bin) file according to the Pagoda2 walkthrough: https://github.com/hms-dbmi/pagoda2/blob/master/vignettes/ pagoda2.walkthrough.oct2018.md. This binary file (Supplementary file 7) can be uploaded to the graphical user interface for exploration of our dataset. We generated a file of cluster labels (for LEC, vEC, VSMC, aEC1, aEC2, aEC3, aEC4, aEC5, aEC6 as well as all other cluster grouping used for analysis), which can also be uploaded to the graphical user interface for visualization of these clusters (Supplementary file 8).

Immunohistochemistry
Mice were euthanized as described above, then perfused through the left ventricle with 5 mL PBS immediately followed by 10 mL ice-fold 4% paraformaldehyde (PFA) in PBS. The left ventricle was then perfused with 6 mL PBS. After aorta dissection, remaining fat tissue was removed with the aorta suspended in PBS in a polystyrene dish. For sectioning, intact thoracic aortae were additionally postfixed in 4% PFA at 4˚for 10 min, then briefly washed with PBS three times. Aortae were then cryoprotected in 30% sucrose in PBS for 2 hr at 4˚C, embedded in a 1:1 mixture of 30% sucrose PBS: OCT over dry ice, then sectioned at 14 mM intervals using a cryostat (Leica Biosystems). E11.5 embryos were fixed in 4% PFA at 4˚for 1 hr, briefly washed with PBS three times, cryoprotected in 30% sucrose in PBS for 24 hr at 4˚C, and embedded as described above for sagittal sectioning with a cryostat.
For staining of mesenteric vessels, mice were euthanized as described above then perfusion-fixed with 1% PFA, followed by post-fixation of mesenteries for 1 hr in 4% PFA. Staining was performed as described for aorta whole-mount preparations with the following exceptions: PBS-T contained 0.2% Triton X-100, blocking solution contained 3% BSA and 5% normal donkey serum in PBS-T, and secondary antibody incubation occurred overnight.

Confocal microscopy and image analysis
Images were acquired using a Zeiss LSM810 confocal microscope equipped with an Plan-Apochromat 20x/0.8 or a Plan-Apochromat 40x/1.4 oil DIC objective. Images were captured using Zen2.1 (Zeiss) software and processed with Fiji (NIH). Zen2.1 software was used to threshold GFP signal and manually count GFP+ nuclei per field ( Figure 6). Fiji was used to quantify GFP signal over PROX1 +ASMA+ and PROX1+ASMA-areas (Figure 8), TSP1 signal normalized to total VEC pixels per field ( Figure 7F

S1P analysis
Plasma S1P was extracted as previously decribed (Frej et al., 2015) with minor modification. Plasma aliquots (5 or 10 mL) were first diluted to 100 mL with TBS Buffer (50 mM Tris-HCl pH 7.5, 0.15 M NaCl). S1P was extracted by adding 100 mL precipitation solution (20 nM D7-S1P in methanol) followed by 30 s of vortexing. Precipitated samples were centrifuged at 18,000 rpm for 5 min and supernatant were transferred to vials for LC-MS/MS analysis (see below).
C18-S1P (Avanti Lipids) was dissolved in methanol to obtain a 1 mM stock solution. Standard samples were prepared by diluting the stock in 4% fatty acid free BSA (Sigma-Aldrich) in TBS to obtain 1 mM and stored at À80˚C. Before analysis, the 1 mM S1P solution was diluted with 4% BSA in TBS to obtain the following concentrations: 0.5 mM, 0.25 mM, 0.125 mM, 0.0625 mM, 0.03125 mM, 0.0156 mM, and 0.0078 mM. S1P in diluted samples (100 mL) were extracted with 100 mL of methanol containing 20 nM of D7-S1P followed by 30 s of vortexing. Precipitated samples were centrifuged at 18,000 rpm for 5 min and the supernatants were transferred to vials for LC-MS/MS analysis. The internal deuterium-labeled standard (D7-S1P, Avanti Lipids) was dissolved in methanol to obtain a 200 nM stock solution and stored at À20˚C. Before analysis, the stock solution was diluted to 20 nM for sample precipitation.
Quantitative linearity was determined by plotting the AUC of the standard samples (C18-S1P) normalized by the AUC of internal standard (D7-S1P); (y) versus the spiked concentration of S1P (x). Correlation coefficient (R2) was calculated as the value of the joint variation between x and y. Linear regression equation was used to determined analyte concentrations.
. Supplementary file 4. scRNA-seq of GFP high and GFP low MAECs. Tabs include transcripts enriched and depleted in each cluster with details at the top of each page. Also included are 2 tabs of transcripts that are commonly enriched or depleted in LEC/vEC/aEC1 (referred to in Figure 3-figure  supplement 4), as well as a tab detailing the intersection of aEC1-enriched transcripts and those differentially expressed in the S1pr1 ECKO dataset ( Figure 6A and B).
. Supplementary file 5. Tabs include Gene Set Enrichment Analysis (GSEA) results using transcripts enriched in each cluster as inputs. Also included is the complete list of sphingolipid-related genes queried (referred to in Figure 3-figure supplement 4A). The list of transcription factors is limited to those identified as aEC1-enriched and -depleted that also had minimal count thresholds after Pagoda2 filtering during hierarchical differential expression analysis.
. Supplementary file 6. Detailed intersections of S1pr1 ECKO up-and down-regulated transcripts according to their cluster assignments from scRNA-seq.
. Transparent reporting form

Data availability
Sequencing data and processed files have been deposited in GEO under the accession GSE139065.
The following dataset was generated:

Author(s)
Year Dataset title Dataset URL