LESR2 is a lymphatic endothelial-specific lncRNA that governs cell proliferation and migration through KLF4 and SEMA3C

Recent studies have revealed the importance of long noncoding RNAs (lncRNAs) as tissue-specific regulators of gene expression. There is ample evidence that distinct types of vasculature undergo tight transcriptional control to preserve their structure, identity, and functions. We determined, for the first time, the global lineage-specific lncRNAome of human dermal blood and lymphatic endothelial cells (BECs and LECs), combining RNA-Seq and CAGE-Seq. A subsequent genome-wide antisense oligonucleotide-knockdown screen of a robust set of BEC- and LEC-specific lncRNAs identified LESR2 as a critical gatekeeper of the global LEC transcriptome. Deep RNA-DNA, RNA-protein, and phenotype rescue analyses revealed that LESR2 acts as a nuclear trans-acting lncRNA modulating, via key epigenetic factors, the expression of essential target genes, including KLF4 and SEMA3C, governing the growth and migratory ability of LECs. Together, our study provides new evidence supporting the intriguing concept that every cell type expresses precise lncRNA signatures to control lineage-specific regulatory programs.


Introduction
The blood and lymphatic vascular systems are essential for the efficient transport of oxygen, nutrients, signaling molecules, and leukocytes to and from peripheral tissues, the removal of waste products, and the preservation of fluid homeostasis. Increased activation or impaired function of these vascular networks represent a hallmark of many pathological conditions, including cancer progression, chronic inflammatory diseases, and diseases leading to blindness [1][2][3] .
During development, the blood vascular system arises from endothelial cell progenitors that differentiate from mesodermal cells, mostly through the expression of the transcription factor (TF) ETV2. Activation of the VEGFA/VEGFR2 signaling and expression of blood vascular endothelial cell (BEC) markers, such as NRP1 and EphrinB2, further differentiate these precursor cells into BECs, which then form the hierarchical network of blood vessels 4 . In contrast, lymphatic vasculogenesis starts after the establishment of the blood circulatory system. Thereafter, a distinct subpopulation of endothelial cells lining the cardinal vein starts differentiating by expressing the TF PROX1, the master regulator of lymphatic endothelial cell (LEC) identity, via the TFs SOX18 and COUPTFII. Once exiting the veins, LECs starts expressing other lymphatic-specific markers, such as podoplanin, VEGFR3, and NRP2, and they migrate, in a VEGFC-dependent manner, to form the primary lymph sacs from which the lymphatic vascular system further develops following sprouting, branching, proliferation, and remodeling processes 5 . However, a nonvenous origin of LECs has also been described in the skin, mesenteries, and heart [6][7][8] . In adulthood, while the blood and lymphatic vasculature are generally quiescent, they can be readily activated in pathological conditions such as wound healing, inflammation, and cancer by disturbance of the natural balance of pro-and anti-(lymph)angiogenic factors 1,9 . Therefore, this complex regulatory network requires precise control of gene expression patterns at both transcriptional and post-transcriptional levels in order to ensure proper maturation, differentiation, and formation of blood and lymphatic vessels.
In this scenario, many studies have recently revealed the importance of a new member of the noncoding RNA clade, termed long noncoding RNAs (lncRNAs), in the regulation of gene activity 10,11 .
In particular, the FANTOM (functional annotation of the mammalian genome) consortium pioneered the discovery of the noncoding RNA world by providing, through Cap analysis of gene expression (CAGE-Seq), the first evidence that large portions of our genome are transcribed, producing a multitude of sense and antisense transcripts 12 . In the latest genome annotation, lncRNAs, which are arbitrarily defined as noncoding RNAs longer than 200 nucleotides, constitute approximately 72% of the transcribed genome 13 , whereas mRNAs comprise only 19%, indicating the need for functional annotation of lncRNAs. Importantly, lncRNAs have recently been shown to display a higher tissuespecificity than mRNAs, suggesting them as new players in the regulation of cell-type-specific gene expression programs 14 .
Since lncRNAs lack a protein-coding role, their primary categorization is based on their genomic location and orientation relative to protein-coding genes 15 . lncRNAs can reside either between protein-coding genes (intergenic, lincRNAs), between two exons of the same gene (intronic lncRNAs), antisense to protein-coding transcripts (antisense lncRNAs), or in promoters and enhancers (natural antisense transcripts or transcribed from bidirectional promoters) [16][17][18] . lncRNAs may regulate gene expression through a multitude of mechanisms depending on their subcellular localization. For instance, in the nucleus, lncRNAs can act as a scaffold for TFs, chromatin remodeling complexes, or ribonucleoprotein complexes (RNPs), indicating a potential role in transcriptional regulation 19 . Nuclear lncRNAs can furthermore act in cis or trans to regulate gene expression by the recruitment of activating and repressive epigenetic modification complexes. Cis-acting lncRNAs, such as the 17-kb X chromosome-specific transcript Xist, regulate gene expression of adjacent genes by directly targeting and tethering protein complexes 20,21 . On the other hand, trans-acting lncRNAs, such as the HOTAIR lncRNA, regulate gene expression at distinct genomic loci across the genome by serving as a scaffold that assists the assembly of unique functional complexes 22 .
In the context of the international FANTOM6 project, which aims to functionally annotate all lncRNAs present in our genome, we first determined lineage-specific lncRNAs associated with human primary dermal LECs and BECs by combining RNA-Seq and CAGE-Seq analyses. Genome-wide functional interrogation after antisense-oligonucleotide (ASO) knockdown of robustly selected LEC and BEC lncRNAs, allowed us to identify a lymphatic endothelial-specific lncRNA, LESR2, to function in the transcriptional regulation of LEC growth and migration. We demonstrated that LESR2 is a trans-acting lncRNA that acts as a protein scaffold in order to facilitate the assembly of unique functional epigenetic complexes involved in gene expression regulation. Through these interactions, LESR2 controls intricated transcriptional networks to fine-tune the expression, above all, of essential proliferation-and migration-related genes, including the tumor-suppressor TF KLF4 and the semaphorin guidance molecule SEMA3C.

Identification of a core subset of vascular lineage-specific lncRNAs
To identify vascular lineage-specific lncRNAs, we performed both RNA-Seq and CAGE-Seq 31 of total RNA isolated from neonatal human primary dermal LECs and BECs (Supplementary Figure 1a, b).
Compared to RNA-Seq, CAGE-Seq allows mapping transcription start sites (TSSs) after quantification of the expression of 5'-capped RNAs 32 . To ensure endothelial cell specificity, we included RNA-Seq and CAGE-Seq data from neonatal human primary dermal fibroblasts (DFs) 33 . In a first step, we performed differential expression (DE) analysis of RNA-Seq and CAGE-Seq of LECs against BECs, LECs against DFs, and BECs against DFs using EdgeR 34 . From defined LEC-or BEC-specific genes (see Methods section), we selected genes annotated as lncRNAs in the recently published FANTOM CAT database 13 . Finally, we overlapped the RNA-Seq and CAGE-Seq results to select lncRNAs identified as differentially expressed using both techniques (Figure 1a The intersection of these two methods revealed 142 LEC and 160 BEC lncRNAs specifically expressed in either LECs or BECs. We defined these subsets as LEC and BEC core lncRNAs ( Figure   1d and Supplementary Table 1).
To characterize the identified LEC and BEC core lncRNA subsets, we analyzed their genomic classification related to protein-coding genes, using the FANTOM CAT 13 annotations. We found that the largest fraction of both LEC and BEC core lncRNAs were categorized as intergenic lncRNAs (47.9% for LEC and 58.8% for BEC), with a significant enrichment compared to all expressed lncRNAs (fold enrichment = 1.9 resp. 2.17, P-value < 0.05) (Figure 1e, f). Gene ontology (GO) analysis of lncRNAs flanking protein-coding genes using Genomic Regions Enrichment of Annotations Tool (GREAT) 35 and g:Profiler 36 showed that both core lncRNA subsets mainly reside near genes related to vascular development, tissue morphogenesis, and endothelial cell function, including proliferation, migration, and adhesion (Supplementary Figure 1c-f). These results are intriguing since several intergenic lncRNAs have previously been reported to play a prominent role in the regulation of gene expression in a cell-specific manner 11 .

Identification of lncRNA candidates for functional characterization by antisense oligonucleotides (ASOs)
To further select lncRNA candidates for genome-wide functional screening, we relied on the FANTOM CAT annotations 13 . Firstly, we filtered for lncRNAs with a conserved transcription initiation region (TIR) and/or exon regions, based on overlap with predefined genomic evolutionary rate profiling (GERP) elements 37 . Secondly, we selected for actively-transcribed lncRNAs with an overlap between TSSs and DNase hypersensitive sites (DHSs). Thirdly, filtering for expression levels in LEC and BEC RNA-Seq and CAGE-Seq datasets (see Methods section) led to the identification of 5 LEC and 12 BEC lncRNAs that are potentially conserved at the sequence level, actively transcribed, and robustly expressed in the respective endothelial cell types (Figure 2a-c). Finally, we identified through qPCR 2 LEC (RP11-536O18.1 and LINC01197), and 2 BEC (LINC00973 and LINC01013) lncRNAs that were consistently differentially expressed between LECs and BECs derived from newborn and adult skin samples (Figure 2d, e). Given their specificity, we renamed these lncRNAs lymphatic endothelialspecific lncRNA 1 and 2 (LESR1 -RP11-536O18.1; LESR2 -LINC01197), and blood vascular endothelial-specific lncRNA 1 and 2 (BESR1 -LINC00973; BESR2 -LINC01013).
Since antisense oligonucleotide (ASO) GapmeRs are more effective in reducing the expression of nuclear lncRNAs than short interference RNAs (siRNAs) 38,39 , we analyzed the subcellular localization of the 2 LEC and 2 BEC lncRNAs in LECs and BECs, using cellular fractionation followed by qPCR.
For LEC lncRNAs, LESR1 was almost equally distributed between cytoplasm and nucleus, whereas LESR2 showed a higher nuclear distribution (Figure 2f). Both BESR1 and BESR2 were mainly localized in the nucleus ( Figure 2g). Therefore, we next used the ASO-based approach to analyze the genome-wide transcriptional changes upon knockdown of the 2 LEC and 2 BEC lncRNAs. After testing their knockdown efficiencies, we selected three out of five ASOs for each lncRNA target  Table 2).

Transcriptional profiling after LESR2-ASOKD indicates potential functions in cell growth, cell cycle progression, and migration of LECs
To investigate the potential functional relevance of the 2 LEC and 2 BEC lncRNAs, we first transfected LECs and BECs with three independent ASOs per target, followed by CAGE-Seq ( Figure   3a and Supplementary Figure 3a, b). Next, we performed DE analysis by comparing the combined results of the ASOs per target with their scramble controls, using EdgeR with a generalized linear model (GLM) 34 . Finally, we defined DE genes by a false discovery rate (FDR) < 0.05 and a log2 fold change (log2FC) > 0.5 resp. < -0.5 (Supplementary Table 3). We found that ASO knockdown (ASOKD) of LESR1 in LECs and of BESR1 and BESR2 in BECs showed rather modest changes in gene expression. LESR1-ASOKD caused changes of only 9 genes (4 up and 5 down), BESR1- genes have previously been reported to play prominent roles in vascular development and differentiation pathways, including PTGS2, KLF4, VEGFA, and ANGPT2 among the upregulated genes, and PROX1, CCBE1, SEMA3C, and ROBO1 among the downregulated genes [40][41][42][43][44] . GO analysis for biological processes using g:ProfileR 36 revealed that indeed both up-and downregulated genes were enriched (P-value < 0.05) for terms related to vascular development. In addition, upregulated genes were mainly involved in cell death, inflammatory signaling, and response to external stimuli, whereas downregulated genes were primarily related to the regulation of cell migration and chemotaxis (Figure 3c and Supplementary Table 4). Gene set enrichment analysis (GSEA) 45 also identified significant (FDR < 0.05) biological processes related to cell migration, chemotaxis, and response to external stimuli/virus. More importantly, several downregulated biological processes were associated with cell growth, cell cycle progression, and cytoskeleton organization (Figure 3d and Supplementary Table 5).
To identify TFs potentially affected by LESR2-ASOKD, we performed motif activity response analysis (MARA) 46 by analyzing the activity of 348 regulatory motifs in TF sites in the proximal promoters of highly expressed genes in knockdown and control samples (see Methods section). We found 19 upregulated and 7 downregulated motifs, among which were binding sites related to several TFs known to be essential for LEC biology, including STAT6, KLF4, NR2F2 (COUPTF-II), and MAFB 43,47,48 (P-value < 0.05, Supplementary Table 6). Interestingly, KLF4 was the only TF to be also upregulated on the transcriptional level upon LESR2-ASOKD. Based on the MARA analysis, we next reconstructed a gene regulatory network with the 255 genes affected by LESR2-ASOKD. We identified modules of up and downregulated genes linked with the identified TF motifs. Among these modules, we found genes associated with endothelial cell proliferation and migration, such as VEGFA, MAFF, ANGPT2, RASD1, PROX1, SEMA3C, and ROBO1 40,49,50 (Figure 3e, f). Overall, these results suggest that the absence of LESR2 had a critical impact on the global transcriptome of LECs by affecting complex TF regulatory networks targeting essential genes largely involved in endothelial cell differentiation, proliferation, and migration.
Before characterizing the biological and molecular function of LESR2 in LECs, we sought to validate its lymphatic endothelial cell-selective expression in situ. Hence, we analyzed the expression levels of LESR2 and specific blood and lymphatic markers in freshly isolated LECs and BECs from human skin samples, using flow cytometry followed by qPCR (Supplementary Figure 3h). As shown in Figure 2d, LESR2 was more highly expressed in LECs isolated ex vivo than BECs. Remarkably, we observed that the LEC specificity of LESR2 was more pronounced in freshly isolated ECs than cultured ECs, similar to the LEC lineage-specific TF PROX1 (Figure 3g).

Knockdown of LESR2 reduces cell growth, cell cycle progression, and migration of LECs in vitro
To investigate the potential effects of LESR2-ASOKD on LEC growth, we performed cell growth assays based on dynamic imaging analysis. We found that LESR2-ASOKD strongly reduced cell growth of LECs over time (Figure 4a and Supplementary Figure 4a, b). To study whether the cell growth phenotype was not due to off-target effects of the ASOs, we also performed cell growth assays after CRISPR interference (CRISPRi) 51 . Consistently, we found that CRISPRi-KD of LESR2 also significantly reduced the growth rate of LECs. However, due to the lower knockdown efficiency, the effect was less prominent compared to ASOKD (Figure 4b and Supplementary Figure 4c Since the transcriptional studies also indicated a potential role of LESR2 in cell migration, we performed wound-closure assays ("scratch assays") after LESR2-ASOKD in LECs. We observed a significant reduction of LEC migration compared to control ( Table 7). Subsequently, we overexpressed the most abundant transcript variant in LECs using a lentiviral vector and analyzed the cell cycle progression and cell migration after LESR2-ASOKD. We performed both assays with the most effective LESR2-ASO2, which binds to the first intron recognizing exclusively the endogenous but not ectopically overexpressed LESR2 LESR2 is a nuclear lncRNA interacting in trans with DNA regions near a subset of differentially expressed genes A first step to study the molecular mechanism of a lncRNA of interest is to analyze its subcellular distribution at the single-molecule level 52 . To this end, we performed single-molecule RNA fluorescence in situ hybridization (smRNA-FISH) in cultured LECs and human skin samples 53 .
Consistent with the cellular fractionation data (Figure 2f), LESR2 was predominantly localized in the nucleus of cultured LECs, showing a broad nuclear distribution with distinct foci (Figure 5a, b). This localization pattern was also observed in lymphatic vessels in human skin (Supplementary Figure 5a), suggesting that LESR2 might exert a chromatin-related function in vitro as well as in vivo.
To further elucidate the possible interactions between LESR2 and chromatin, we performed chromatin isolation by RNA purification followed by DNA sequencing (ChIRP-Seq) 54 . Cross-linked LECs were hybridized with two biotinylated probe sets (Odd and Even, internal control) tiling LESR2 (Supplementary Table 8). Probes targeting LacZ were used as an additional control. After pull-down, the percentage of retrieved RNA was assessed (Supplementary Figure 5b), and DNA was subjected to sequencing. Using a previously published analysis pipeline 54 , we found 2,258 binding sites of LESR2 to be at least 3-fold significantly enriched compared to input (P-value < 0.05; see Methods section), including a peak in the LESR2 exon one region as pull-down control (Supplementary Figure   5c and Supplementary Table 9).
To identify candidate genes directly regulated by LESR2, we first analyzed the genomic distribution of LESR2 binding sites. Out of 2,258 binding sites, 1,497 mapped within protein-coding genes (65.5%), with a large fraction residing in introns (1,010 peaks, 68.3%) (Figure 5c). Since only 19% of all annotated genes are categorized as protein-coding in the FANTOM CAT database 13 , these results suggested a preference of LESR2 to interact with regulatory regions near protein-coding genes (fold enrichment = 3.24, P-value < 0.05). Therefore, we focused on the identified 1,607 protein-coding genes displaying at least one LESR2 binding site within their promoters, exons, or introns. LECs for downregulated genes, or in BECs for upregulated genes, implicating these genes as potential downstream targets of LESR2 (Figure 5f).

KLF4 and SEMA3C
Among the 44 potential downstream targets of LESR2, KLF4 caught our attention as potential cell proliferation regulator given its well-established tumor suppressor role 55 and the previously observed upregulation at the RNA level as well as increased TF binding activity upon LESR2 knockdown ( Figure 3). Among cell migration regulatory molecules, we focused on one member of the semaphorin protein family, SEMA3C, that was previously shown to enhance migration in endothelial cells 56 .
To functionally characterize the relationship between LESR2 and KLF4 as well as SEMA3C, we performed the experimental strategies represented in Figure 6a. For KLF4, we analyzed the cell cycle progression of LECs after LESR2-ASOKD, followed by siRNA knockdown of KLF4. As expected, LESR2-ASOKD resulted in an upregulation of KLF4 as well as an increase of G0-arrested LECs.

LESR2 interacts with several protein complexes to exert its regulatory function on gene expression
To identify proteins that are potential co-regulator of LESR2 target genes, we performed in vitro biotin-LESR2 pull-down assays 57 . Nuclear extracts of LECs were incubated with the biotinylated full-length LESR2 transcript and its antisense as negative control (Supplementary Figure 7a, b). After streptavidin bead separation, mass spectrometry was performed to identify possible interacting proteins. Initial analysis identified a total of 642 proteins. After filtering for proteins present in both replicates but absent in the antisense control, we found 59 proteins to interact with LESR2 (Supplementary Table 10). GO analysis for molecular functions and cellular compartments using g:ProfileR 36 confirmed that the 59 identified proteins were significantly enriched for nuclear RNAbinding proteins (Supplementary Figure 7c, d). Protein-protein interaction analysis using STRING 58 revealed that a large fraction of these 59 proteins were associated with RNA-processing functions, such as RNA splicing, RNA polyadenylation, and RNA nuclear transport. Furthermore, six proteins were associated with chromatin remodeling and three with nuclear organization, suggesting that LESR2 may operate at several levels to regulate gene expression (Figure 7a).
To screen for protein candidates, we analyzed the RNA expression of the 59 proteins interacting with LESR2 in LECs versus BECs. Four proteins (DDX39A, NUMA1, RBBP7, and DDX5) had a logFC greater than 0.5 in LECs and a unique peptide detection greater than 5 ( Figure 7b). Among these proteins, we identified the histone-binding protein RBBP7, which has previously been reported to be involved in the regulation of many cellular functions, including proliferation and migration 59,60 .
Subsequent RNA immunoprecipitation assays in LECs validated the interaction between RBBP7 and LESR2, suggesting RBBP7 as a potential mediator of LESR2 gene regulatory functions (Figure 7c, d).
In summary, our multilayered mode of action analysis demonstrates that LESR2 is a nuclear lncRNA that interacts with essential epigenetic partners to regulate cell growth and cell migration of LECs by tuning the expression of distinct target genes, in particular KLF4 and SEMA3C (Figure 7e).

Discussion
Precise regulation of proliferation, migration, and maintenance of cellular identity are not only essential to ensure proper development and integrity of the vascular systems, but also to guarantee that LECs and BECs are able to perform their necessary functions 5 . In this study, we characterized, for the first time, the global lineage-specific lncRNAome of LECs and BECs, and analyzed the transcriptional impacts after ASO-mediated knockdown of LEC-and BEC-specific lncRNAs followed by CAGE-Seq. Importantly, we identified LESR2 as a novel lymphatic-specific lncRNA that is essential in the regulation of lymphatic endothelial cell growth and migration.
By integrating RNA-Seq and CAGE-Seq transcriptome profiling, we showed that LECs and BECs express a specific cohort of lncRNAs, mainly residing near vascular-related protein-coding genes.
These results are in accordance with the intriguing concept that cells might display a set of lncRNAs explicitly expressed to function in the fine-tuning of cell-type-specific gene expression programs 11 .
Most notably, our selection strategy highlighted 2 LEC and 2 BEC lncRNAs that are robustly and differentially expressed in the respective endothelial lineage. These candidates therefore represent the first set of lineage-specific lymphatic and blood vascular endothelial cell lncRNA markers.
The nuclear localization of our lncRNA candidates coincides, to some extent, with previous findings, demonstrated by RNA in situ hybridization, that lncRNAs are commonly located in the nucleus 61 . To investigate the biological functions of these lncRNA candidates, we used the antisense oligonucleotide GapmeR knockdown approach, given its higher efficiency in targeting nuclear RNA transcripts over siRNA 38 . Additionally, our ASO design strategy proved to be very successful, resulting in an extremely high knockdown efficiency of all four targets. In the general experimental design, we have not taken into consideration the use of CRISPRi due to several practical limitations to study lncRNAs. For instance, as opposed to ASOKD, CRISPRi may also interfere with the transcription of overlapping or neighboring transcriptional units, and it is not able to distinguish the cis and transacting functions of lncRNAs 62,63 . However, we are aware that CRISPRi provides less pronounced offtarget effects compared to other loss-of-function methods 64 .
Importantly, our study identified the first LEC lncRNA with specific biological functions. Through a multilayered analysis, including differential expression analysis, GO, GSEA, and MARA, we found that LESR2, originally annotated as LINC01197, is a critical gatekeeper of the global transcriptome of LECs by influencing complex TF regulatory networks regulating essential targets largely involved in the control of LEC growth and migration. These "molecular phenotypes" observed after LESR2 knockdown were confirmed in vitro by analyzing the cell growth profile, cell cycle progression, and wound closure ability of primary human LECs. As shown in a previous study 33 , we thus were able to distinguish LESR2 as a bona fide functional lncRNA in LECs by combining sequencing data analyses and cellular phenotype assays.
The nuclear localization of LESR2, as shown by subcellular fractionation and smRNA-FISH, hinted a chromatin-related function. Indeed, as confirmed by RNA-DNA interaction assay, we revealed that LESR2 interacts, predominantly in trans, with DNA regions near a subset of differentially expressed genes. In addition, RNA-protein interaction assays indicated a potential scaffold function of LESR2 in recruiting proteins involved in several levels of gene expression regulation, including chromatin organization. These results are in line with the general model in which lncRNAs are crucial for the assembly of unique protein complexes and for guiding them to specific target sites 65 . Specifically, we found that LESR2 interacts with RBBP7, a protein previously reported to be part of several multiprotein complexes that are involved in chromatin remodeling, histone post-transcriptional modification, and gene expression regulation [66][67][68] . Intriguingly, RBBP7 is a relevant constituent of the polycomb repressive complex 2 (PRC2) complex, which was also previously shown to interact with 20% of the lncRNAs in human cells 69 . It is conceivable that LESR2 might act as an epigenetic regulator to recruit or guide protein partners to influence the three-dimensional structure of the genome. In this setting, the differential CAGE-Seq peak intensities at the target TSSs after LESR2-ASOKD might provide a hint on the function of LESR2 in mediating the transcriptional machinery access at the site of transcription. In fact, significant correlations at TSS regions between RNA polymerase II occupancy and CAGE-Seq signal have been reported 70,71 .
We showed that LESR2 exerts its effects on LEC functions, at least in part, via the modulation of KLF4 and SEMA3C, since knockdown of KLF4 or overexpression of SEMA3C partially restored the cellular phenotypes observed upon LESR2-ASOKD. Previous studies have reported that KLF4 is a tumor-suppressor TF that, once upregulated, inhibits cell growth and induces cell cycle arrest [72][73][74] . A primary mechanism by which KLF4 regulates cell growth is via the induction of CDKN1A expression, a gene encoding a cyclin-dependent kinase (CDK) inhibitor 73,74 . Consistently, we found that CDKN1A was also upregulated upon knockdown of LESR2 (Supplementary Table 3), suggesting that suppression of the KLF4-CDK1NA axis through LESR2 is required for the maintenance of a proliferative state of LECs. Moreover, MARA analysis highlighted that the trans-acting activity of LESR2 has a general inhibitory effect on the binding activity of KLF4, suggesting an additional genome-wide interplay between KLF4 and LESR2 to modulate sensitive targets indispensable for LEC function.
A recent study reported that viral ectopic expression of FOXO1 in ductal pancreatic adenocarcinoma cells (PDAC) upregulates LINC01197 (LESR2), resulting in the inhibition of cell proliferation 75 .
Interestingly, FOXO1 has previously been shown to promote LEC migration and to participate in the regulation of lymphatic development 76,77 . Moreover, in our sequencing data, FOXO1 was significantly more highly expressed in LECs than BECs.
SEMA3C belongs to the semaphorin class 3 guidance cue molecules, which mainly bind to a receptor complex composed of neuropilins (NRP1 or NRP2) and plexins (PLXNA1-A4 and PLXND1) 78 . Our findings that overexpression of SEMA3C partially rescued the LESR2 inhibition of LEC migration are congruent with previous reports showing that SEMA3C has pro-migratory activities in several cell types [79][80][81] , including endothelial cells 82 . In support of this claim, our sequencing data revealed that LECs expressed the two neuropilins, NRP1 (at a low level) and NRP2, as well as the plexins A1-A4 and D1. Additionally, knockdown of LESR2 also significantly reduced the expression of plexin A4 (Supplementary Table 3), pinpointing LESR2 as an intermediate player in semaphorin signaling.
While we initially identified LESR2 as a LEC-specific lncRNA by sequencing of cultured human LECs, smRNA-FISH and FACS validated its lymphatic specificity in human skin in situ. Future studies are needed to investigate its expression pattern and mechanistic role in pathological conditions associated with impaired lymphatic function (e.g., lymphedema), or active lymphangiogenesis (e.g., tumor metastasis and wound healing).
It is of interest that the knockdown of three other lncRNA candidates showed merely a minor or no impact on the transcriptome of LECs and BECs after ASO-mediated knockdown. Likely, these could be due to four potential reasons. First, they may have alternative functions unrelated to transcriptional regulation, such as ribozymes or riboswitches 83 and translation initiation regulators 84 . Second, the act of transcription, rather than the lncRNA product of this transcription, may be functional by having, for instance, an enhancer-like function 85 . Third, they may function as molecular signals at a specific time and place in response, for example, to unique stimuli 86 . Finally, although clearly differentially expressed, all three lncRNA candidates might not be functional and might just be part of transcriptional noise 87 . Therefore, future research is needed to elucidate the biological role and function of these lncRNAs in LECs or BECs.

Taken together, our study enumerates the collection of lncRNAs explicitly expressed in LECs and
BECs and highlights, through the functional characterization of LESR2, the importance of those lncRNAs in the regulation of lineage-specific endothelial cell functions.

Isolation of adult primary skin LECs and BECs from biopsies
LECs and BECs were obtained from the abdominal or breast skin of healthy adult subjects admitted for plastic surgery at the University Hospital Zurich. Written informed consent was obtained from each donor/tissue collection, as approved by the Ethics Committee of the Kanton Zurich (2017-00687).

Cell culture
Primary human dermal LECs and BECs were isolated from neonatal human foreskin as described previously 88 and cultured in endothelial basal medium (EBM, Lonza) supplemented with 20% FBS, 100U/mL penicillin and 100µg/mL streptomycin (Pen-Strep, Gibco), 2mM L-glutamine (Gibco), 10µg/mL hydrocortisone (Sigma) on sterile dishes/plates (TPP) pre-coated with 50µg/mL purecol type I bovine collagen solution (Advanced BioMatrix) in DPBS at 37°C in a 5% CO2 incubator. LECs were additionally cultured in the presence of 25µg/mL cAMP (Sigma); BECs in the presence of endothelial cell growth supplement ECGS/H (PromoCell). Cells were used at passage 6-7 for RNA-Seq and CAGE-Seq experiments, and at passage 8-9 in in vitro and biochemistry experiments. Primary human dermal LECs and BECs were isolated from adult human skin as described above and cultured in EGM-2 complete medium in vessels coated with fibronectin diluted in DPBS and cells were cultured at 37°C in a 5% CO2 incubator. Cells were used between passage 2 and 6. HEK293T cells were cultured in DMEM with glutamax (Gibco) supplemented with 10% FBS (Gibco) and Pen-Strep at 37°C in a 5% CO2 incubator. All cells were routinely tested for mycoplasma contamination using the MycoScope PCR Mycoplasma Detection kit (Genlantis).

RNA isolation, reverse transcription, and qPCR
If not differently specified, total RNA was isolated using the NucleoSpin RNA kit (Machery Nagel),  Table 11.

Western blot analysis
To perform western blot analysis, the protein concentration of lysates was first determined using the Microplate BCA protein assay kit -reducing agent compatible (Thermo Fisher Scientific), according to the manufacturer's instruction. To 5-30µg of total protein, SDS sample buffer and reducing agent (Thermo Fisher Scientific) were added to a 1x final concentration. Then, samples were heated for 5min at 70°C and size-separated by electrophoresis using 1.5mm 4-12% NuPAGE Bis-Tris protein gels and 1x NuPAGE MES SDS running buffer (Thermo Fisher Scientific) for 35-50min at 200V.
Proteins were transferred to a nitrocellulose membrane (Merck Millipore) for 1h at 20V. Protein loading was checked by staining membranes with Ponceau staining solution (Sigma) for 2min at room temperature. Membranes were blocked with 5% milk powder in TBST (50mM Trizma Base, 150mM NaCl, 0.1% Tween 20, pH 8.4) for 1.5h at room temperature. Membranes were then stained overnight at 4°C with primary antibodies (see below) diluted in TBST. Blots were washed three times with TBST for 15min at room temperature and subsequently incubated for 2h at room temperature with HRPconjugated secondary antibodies (goat anti-rabbit, Dako; rabbit anti-goat, R&D Systems) at a dilution of 1:5000 in TBST. After washing five times with TBST for 15min at room temperature, blots were developed with clarity western ECL substrate (Bio-Rad) and imaged on a ChemiDoc imaging system (Bio-Rad). All antibodies are listed in the Supplementary Table 11.

Lentivirus production
For production of lentiviruses, 2.5x10 6 HEK293T cells were seeded into 10cm dishes and cultured overnight. One hour before transfection, the medium was replaced with an antibiotic-free medium containing 25µM chloroquine (Sigma). The transfection mixture was subsequently prepared as follows. In a first tube, 1.3pmol psPAX2 (12260, Addgene), 0.72pmol pMD2.G (12259, Addgene), and 1.64pmol of target vector were mixed in 500µL Opti-MEM (Gibco). In another tube, polyethylenimine (PEI, Sigma) was added to 500µL Opti-MEM in a 1:3 ratio to total DNA content. PEI-containing Opti-MEM was transferred dropwise to the plasmid-containing Opti-MEM, and the mixture was incubated for 20min at room temperature. Finally, the transfection mixture was transferred dropwise to the HEK293T cells. 24h post-transfection, the medium was changed with 8mL complete medium containing 10µM forskolin (Sigma). Lentiviruses were harvested after 48h, filtered with a 0.45µm PES filter (TPP), and stored at -80°C. The titer of each virus was determined using the Lenti-X Go-Stix Plus kit (Takara Bio), according to the manufacturer's instruction.

RNA-Seq and CAGE-Seq of primary LECs and BECs
2.5x10 5 LECs and BECs were seeded in duplicates into 10cm dishes and cultured for three days until 70% confluence was reached. At this point, 8mL EBM consensus medium (EBM supplemented with 20% FBS, L-glutamine, and Pen-Strep) was added to both cell types. After 48h, total RNA was harvested and isolated using RNeasy mini kit (Qiagen). DNA digestion was performed using the RNase-Free DNase set (Qiagen). The identity of BECs and LECs was checked by qPCR ( Supplementary Figure 1a, b). Primers are listed in the Supplementary Table 11. LEC and BEC total RNA were then subjected to ribosomal-RNA depleted RNA sequencing (RNA-Seq) and nAnT-iCAGE sequencing (CAGE-Seq) protocols as previously described 31 .

Differential gene expression analysis of CAGE-Seq and RNA-Seq data
For both RNA-Seq and CAGE-Seq data, read alignment was performed, and expression tables were generated as described previously 33 . Next, we performed differential expression (DE) analysis of RNA-Seq and CAGE-Seq of LECs against BECs, LECs against DFs, and BECs against DFs using EdgeR (ver. 3.12.1) 34,89,90 . LEC-associated genes: false discovery rate (FDR) < 0.01 and log2 fold change (FC) LECs/BECs > 1 and log2FC LECs/DFs > 1. BEC-associated genes: FDR < 0.01 and log2FC LECs/BECs < -1 and log2FC BECs/DFs > 1. LEC-and BEC-associated lncRNAs were selected according to their annotation as "lncRNA" in the FANTOM CAT database 13 . LEC and BEC core lncRNAs were then defined as the overlap between RNA-Seq and CAGE-Seq analyses.

Genomic classification and flanking gene analyses of LEC and BEC core lncRNAs
To analyze the genomic classification of LEC and BEC core lncRNAs, we first determined which lncRNAs were broadly expressed in either BECs or LECs by applying an expression level cut-off of 0.5 on both TPM (RNA-Seq) and CPM (CAGE-Seq). Next, LEC/BEC core lncRNAs and broadly expressed lncRNAs were classified according to the "CAT gene class" category of the FANTOM CAT database 13 . Finally, enrichment analysis of LEC/BEC intergenic lncRNAs versus broadly expressed intergenic lncRNAs was performed using SuperExactTest (ver. 1.0.0) 91 . As background, total annotated lncRNA in FANTOM CAT database 13 (n = 90,166) were used.
To determine flanking protein-coding genes, we uploaded lists containing transcriptional start sites (TSSs) of LEC/BEC core lncRNAs to the GREAT webtool (http://great.stanford.edu/public/html/, ver. 4.0.4) 35 . As association rule, we used the "two nearest genes" option with a maximal extension from the lncRNA TSS of 10 Mb. Gene Ontology (GO) analysis was then performed on determined proteincoding genes, using g:Profiler (ver 0.6.7) 36 . Genes with TPM (RNA-Seq) and CPM (CAGE-Seq) < 0.5 were used as background. The gprofiler database Ensembl 90, Ensembl Genomes 37 (rev 1741, build date 2017-10-19) were used. Only GO terms with P-value < 0.05 were used for further analysis.

Antisense oligonucleotide design and efficiency test
Five locked nucleic acids (LNA) phosphorothioate GapmeRs per lncRNA target were designed as described previously 33   Opti-MEM according to the manufacturer's instructions. After 5min incubation at RT, the Lipofectamine-ASO mixture was added dropwise to the cells. LECs and BECs were harvested after 48h post-transfection. For these samples, total RNA was isolated using the RNeasy mini kit. DNA digestion was performed using the RNase-Free DNase set. Knockdown efficiency for each ASO was checked by qPCR, as described above. Samples with at least 50% knockdown efficiency were subjected to CAGE-Seq. Knockdown efficiency was also confirmed after CAGE-Seq (Supplementary

Differential gene expression analysis of ASOKD data
All ASOs for each targeted lncRNA were compared against scramble control ASO (NC_A) libraries from corresponding cell types. Genes with expression >= 5 CPM in at least two CAGE libraries (targeted lncRNA ASOs + scramble control ASO (NC_A) CAGE libraries) were tested for differential expression (DE). A generalized linear model (GLM) was implemented for each targeted lncRNA to perform DE analysis using EdgeR (ver. 3.12.1) 34,89,90 . Genes with |log2FC| > 0.5 and FDR < 0.05 were defined as differentially expressed genes and used for the downstream analysis.

Gene ontology (GO) analysis of ASOKD data
GO analysis was performed separately on upregulated and downregulated genes, using g:Profiler (ver 0.6.7) 36 , as described above. Genes with expression >= 5 CPM in at least two CAGE libraries were used as background. All the significant GO terms (P-value < 0.05) were used for further analysis and are listed in Supplementary Table 4.
Enriched GO biological processes were selected and organized in a network using Cytoscape (ver. 3.6.1) and the plugin Enrichment Map 93 . Gene-set filtering was set as following: FDR q-value cutoff < 0.05 and P-value cut off < 0.001. Gene-set similarity cutoff was set as < 0.5 with an Overlap Metric.
Genes were filtered by expression. Terms were then organized manually according to their biological meaning using the Cytoscape plugin Wordcloud 93 . Significant GO biological processes after GSEA are listed in Supplementary Table 5.

Motif activity response analysis (MARA) of ASOKD data
MARA was performed for BEC and LEC separately using promoter expression for all the knockdown (KD) and scramble control ASO (NC_A) libraries. All promoters with expression >=1TPM at least in 5 CAGE libraries were used for the analysis. MARA was performed as described previously 33 for the motifs in SwissRegulon (released on 13 July 2015) 94 . Student's t-test was performed to identify differentially active motifs due to the lncRNA-ASOKD. Motifs with significant motif activity differences (P-value < 0.05) compared to the controls (NC_A) were used for downstream analysis.
Significant TF motifs after MARA are listed in Supplementary Table 6.
Cloning sgRNA targeting LESR2 and establishment of dCas9 expressing LECs sgRNAs targeting LESR2 were designed using the online CRISPR design tool from the Zhang lab, MIT (http://crispr.mit.edu/). 250bp upstream of the highest CAGE-Seq peak were used as the design region. We then selected then 3 sgRNAs to be cloned into lentiGuide-Puro (52963, Addgene) as previously described 95 . Briefly, each pair of oligos was first annealed and phosphorylated using T4 PNK (New England BioLabs) using the following program: 30min at 37°C, 5min at 95°C, and then  Table 11. Lentiviruses containing pHAGE EF1a dCas9-KRAB (50919, Addgene, with custom blasticidin cassette), scramble control sgRNA, or each of the sgRNA targeting LESR2 were produced as described above.
To establish dCas9-overexpressing LECs, 1.2x10 5 LECs were seeded into pre-coated 6-well plates (TPP) and infected with medium containing dCas9-KRAB lentiviruses diluted at a 10 multiplicity of infection (MOI) and 5µg/mL polybrene (hexadimethrine bromide, Sigma). Plates were then sealed with parafilm and centrifuged at 1,250rpm for 1.5h at room temperature. The next day, the medium was changed, and positively infected cells were selected with 10μg/mL blasticidin (InvivoGen). Once confluent, at least 5x10 5 dCas9-KRAB-expressing LECs were split into pre-coated 10cm dishes and cultured under antibiotic selection until confluency. After checking RNA and protein levels of dCas9-KRAB as described previously 95 , LECs were then used in the cell growth profile experiment.

Cell growth profiling after LESR2-ASOKD and -CRISPRi
For cell growth profiling after ASOKD, 3,000 LECs per well were seeded into a 96-well plate and cultured overnight. LECs were then transfected with 20nM of scramble control ASO or three ASOs targeting LESR2 (Exiqon) and 0.2µL Lipofectamine RNAiMAX previously mixed in 20µL Opti-MEM according to the manufacturer's instructions.
For cell growth profiling after CRISPRi, 3,000 dCas9-expressing LECs per well were seeded into a 96-well plate and grown overnight. LECs were then infected with 50 MOI of lentiviruses containing vectors expressing scramble control sgRNA or three sgRNAs targeting LESR2 diluted in complete growth EBM medium supplemented with 5µg/mL polybrene. After 24h, the virus-containing medium was changed.
In both experiments, LECs were continuously imaged every 3h over three days with 4 fields per well using the IncuCyte ZOOM live-cell imaging system (Essen Bioscience). Confluence in each well was determined using IncuCyte ZOOM software. The normalized growth rate was calculated as the slope of linear regression and normalized to control. To check knockdown efficiency, LECs were harvested after 72h post-transfection. Total RNA was isolated using the RNeasy mini kit. qPCR was performed using One Step SYBR PrimeScript RT-PCR kit (Takara Bio) on a 7900HT real-time system (Applied Biosystems).

Cell cycle analysis after LESR2-ASOKD
Cell cycle analysis was adapted from 96 . 5x10 5 LECs were seeded into 10cm dishes and cultured overnight in starvation medium (EBM supplemented with 1% FBS). The next day, LECs were transfected with scramble ASO or three ASOs targeting LESR2, as described above. As additional controls, LECs were incubated in the starvation medium as non-proliferative control or treated with 100ng/mL nocodazole (Sigma) as mitosis control. 24h post-transfection, LECs were detached and collected. Floating LECs were also collected in the same tube. 2x10 5 LECs per replicate were then transferred into a 96 U-bottom plate (Greiner bio-one). Aliquots of approx. 1x10 5 LECs were lysed and subjected to RNA isolation, reverse transcription, and qPCR, as described above. After a DPBS wash, LECs were stained with Zombie NIR (BioLegend) diluted 1:500 in DPBS for 15min at room temperature in the dark. Subsequently, LECs were washed with DPBS and fixed in 70% ethanol overnight at -20°C. After permeabilization, LECs were washed twice in DPBS and stained with mouse anti-human Ki-67 antibody (Dako) diluted 1:800 in FACS buffer (DPBS with 1mM EDTA and 2% FBS) for 30min at room temperature in the dark. As a negative control, one sample was stained with mouse IgG isotype control (clone 11711, R&D Systems) diluted 1:5 in FACS buffer. After two washes in DPBS, LECs were then stained with donkey Alexa488-conjugated anti-mouse secondary antibody (Thermo Fisher Scientific) diluted 1:500 in FACS buffer for 30min at room temperature in the dark.
Next, LECs were washed once in DPBS and incubated with 20µg/mL RNase A (Machery-Nagel) diluted in DPBS for 60min at room temperature. Finally, LECs were incubated with 10µg/mL propidium-Iodide diluted in staining buffer (100mM Tris, 150mM NaCl, 1mM CaCl2, 0.5mM MgCl2, and 0.1% Nonidet P-40) for 20min at RT. Flow cytometry was performed with a CytoFlex S instrument (Beckman Coulter). Data were analyzed with FlowJo (ver. 10.lr3). Gating was done using live/dead to gate living cells, isotype control to gate proliferative/non-proliferative cells, starvation medium-treated sample to gate G0 population, and nocodazole-treated sample to gate mitotic cells.

Apoptosis assay after LESR2-ASOKD
3,000 LECs per well were seeded into a 96-well plate and cultured overnight. LECs were then transfected with 20nM of scramble control ASO or three ASOs targeting LESR2, as described above.
As a positive control, additional cells were treated for 2h with 4µM staurosporine (Sigma). 48h posttransfection, LECs were fixed in 3.7% PFA (AppliChem) in DPBS for 20min at room temperature.
After three washes with DPBS, LECs were blocked in blocking buffer (DPBS with 5% donkey serum (Sigma) and 0.3% triton X-100) for 1h at room temperature. LECs were then stained with rabbit anticleaved caspase 3 (Asp175) antibody (Cell Signaling) diluted 1:400 in antibody buffer (DPBS with 1% BSA and 0.3% triton X-100) overnight at 4°C. The next day, LECs were washed three times for 5min in DPBS and subsequently stained with donkey Alexa488-conjugated anti-rabbit secondary antibody (Thermo Fisher Scientific) diluted 1:1000, and Hoechst dye (Thermo Fisher Scientific) diluted 1:2000 in antibody buffer for 1h at room temperature in the dark. After three washes for 5min in DPBS, the plate was imaged using a fluorescence microscope (Zeiss Axiovert 200M), and four images at a 5X magnification were taken for each well. Percentage of cleaved caspase 3-positive cells was determined using a self-built macro developed with ImageJ (ver. 2.0.0-rc-69/1.52i) 97 . To check knockdown efficiency, an extra plate was lysed at the end of the experiment and subjected to qPCR, as described above.
Wound closure assay after LESR2-ASOKD 2.5x10 5 LECs were seeded into 6cm dishes and cultured for at least 6h. LECs were then transfected with 20nM of scramble control ASO or three ASOs targeting LESR2 and 8µL Lipofectamine RNAiMAX previously mixed in 800µL Opti-MEM according to the manufacturer's instructions. 24h posttransfection, LECs were detached, and 15,000 were seeded into 96-well plates and cultured overnight. The next day, LECs were incubated for 2h with 2µg/mL mitomycin C (Sigma) in complete EBM medium. After 2h incubation with proliferation inhibitor, LECs were scratched using a wounding pin replicator (V&P Scientific), according to the manufacturer's instructions, complete EBM medium was replaced, and LECs were incubated for 9h. Images of the scratched areas were taken at 5x magnification and at time points 0h, 3h, 6h, and 9h using a bright field microscope (Zeiss Axiovert 200M). Images were analyzed using TScratch (ver. 1.0) 98 . To check knockdown efficiency, LECs were harvested after 9h, and total RNA, cDNA synthesis, and qPCR were performed as described above.

Trans-well migration assay after LESR2-ASOKD
2.5x10 5 LECs were seeded into 6cm dishes and cultured for at least 6h. Transfection of scramble control ASO or three ASOs targeting LESR2 was performed as described above. 24h posttransfection, 50,000 LECs were seeded in starvation medium into a trans-well (24-well plate, 6.5mm insert, 8µm polycarbonate membrane, Costar) previously coated with collagen-I on both sides of the membrane. To check knockdown efficiency, an aliquot of each condition was lysed and subjected to qPCR, as described above. After 4h incubation at 37°C, LECs were fixed with 3.7% PFA in DPBS for 10min. After a DPBS wash, nuclei were stained with Hoechst dye diluted 1:1000 in DPBS for 10min at room temperature in the dark. LECs on the upper side of the membrane were then removed with cotton swabs, and the membrane was thoroughly washed with DPBS. Finally, membranes were cut and mounted on microscope slides with mowiol 4-88 (Calbiochem). For each membrane, four images were taken at 5x magnification using a fluorescence microscope (Zeiss Axioskop 2 mot plus). Cell migration was quantified by counting the nuclei per image field using a self-built macro developed with ImageJ (ver. 2.0.0-rc-69/1.52i) 97 . GENEzol reagent (Geneaid Biotech). After 5min incubation at room temperature, 200µL chloroform were added to the homogenized nuclear fraction. After vortexing for 10s, the nuclear fraction was spun down at 13,000rpm for 15min at 4°C. The upper aqueous phase was transferred into a new tube. For cytoplasmic RNA, on the other hand, two volumes of phenol:chlorofrom:isoamyl alcohol mixture (PCA, 25:24:1, Sigma) were added to the cytoplasmic fraction and vortexed for 1min. After spinning for 5min at 4,000rpm, the upper aqueous phase was transferred into a tube. One volume of isopropanol was added to both the nuclear and cytoplasmic aqueous phases and mixed by inverting the tubes several times. After 10min incubation at room temperature, the tubes were centrifuged at 13,000rpm for 10min at 4°C. RNA pellets were subsequently washed with 75% ethanol and dried for 10min at room temperature. Dried RNA pellets were resuspended in RNase-free water and incubated for 15min at 58°C on a heating block. Finally, both samples were subjected to cDNA synthesis, and qPCR was performed as described above.

Identification of transcripts variants of LESR2 in LECs
The transcriptional start site (TSS) of LESR2 was determined by examining the CAGE-Seq signal The primer sequence was 5'-GATTACGCCAAGCTTTTGTGAGCCACTGCGTTCT-3', and the annealing temperature was 62°C. Polyadenylation of LESR2 transcripts was determined using qPCR and comparing the expression values between cDNA synthesis with random and oligodT20 primers.
Expression levels of the identified primers were determined using qPCR as described above. Primers to amplify LESR2 transcripts are listed in the Supplementary Table 11. Sequences of LESR2 transcript variants are listed in Supplementary Table 7. 2.0.0-rc-69/1.52i) 97 was used to quantify the nuclear versus cytoplasmic localization of LESR2 by applying a max intensity projection. After determining the nuclear surface using the Hoechst dye signal, plugin "analyze particles" was used to count spots present either in the nuclear or cytoplasmic area.

Simultaneous smRNA-FISH and immunostaining in human skin tissues
Normal human skin samples were obtained from plastic surgery. Immediately after dissection, samples were fixed in fresh 10% neutral buffered formalin (NBF) for 24h at room temperature. Fixed samples were then dehydrated using a standard ethanol series followed by xylene, and were embedded in paraffin. Using a microtome, 5µm skin sections were cut and mounted on superfrost plus slides (Thermo Fisher Scientific). Slides were dried overnight at room temperature. Simultaneous smRNA-FISH and immunostaining was performed using the RNAScope Multiplex Fluorescent Reagent kit v2 (Advanced Cell Diagnostics), according to the manufacturer's instructions and technical note 323100-TN. Briefly, slides were backed for 1h at 60°C and deparaffinized with a series of xylene and ethanol washes. Once dried for 5min at 60°C, slides were incubated with RNAscope hydrogen peroxide for 10min at room temperature. Target retrieval was then performed for 10min in RNAscope target retrieval solution using a steamer constantly held at 94-95°C. A hydrophobic barrier was drawn on each slide and let dry overnight at room temperature. The next day, protease treatment was performed for 30min at 40°C in a HybEZ Oven (Advanced Cell Diagnostics). Afterward, slides were incubated for 2h at 40°C with the following RNA scope probes: Hs-LESR2-C1 (563721-C1, for 30min at room temperature. All slides were finally counterstained with Hoechst dye diluted 1:10000 in DPBS for 5min at room temperature and mounted with proLong gold antifade mountant (Thermo Fisher Scientific). Z-stacks of fluorescence images spanning over the entire tissue section were acquired using an inverted confocal microscope (Zeiss LSM 880). Z-projection of acquired images was done using ImageJ (ver. 2.0.0-rc-69/1.52i) 97 .

FACS sorting of primary BECs and LECs followed by qPCR
Single-cell suspensions from skin samples of adult subjects were prepared as described above.
Subsequently, isolated single cells were stained with mouse anti-human CD34 biotinylated antibody (clone 581, Thermo Fisher Scientific) diluted 5µL for 10 6  After washing in FACS buffer, isolated single cells were filtered and sorted for living, CD45-CD31+ CD34low (LEC), or CD34 high (BEC) directly into test tubes containing 250µL RLT plus lysis buffer, using a FACSAria (BD Biosciences). RNA was isolated using the RNeasy Plus Micro kit (Qiagen), according to the manufacturer's instructions, including DNase digestion. qPCR was performed on an HT7900 system and analyzed as described above.

Chromatin isolation by RNA purification followed by sequencing (ChIRP-Seq)
Chromatin isolation by RNA purification followed by DNA sequencing (ChIRP-Seq) was performed as previously described 54  (Thermo Fisher Scientific) for 30min at 37°C with rotation. After five washes with washing buffer (2xSSC, 0.5% SDS, 1mM PMSF), 100µL beads were used for RNA isolation and 900µL for DNA isolation. RNA aliquots were incubated in 80µL RNA proteinase K solution (100mM NaCl, 10mM Tris-Cl pH 7.0, 1mM EDTA, 0.5% SDS, 0.2U/µL proteinase K (Ambion)) for 45min at 50°C with rotation and boiled for 10min at 95°C. 500µL Qiazol (Qiagen) were added to each sample, and RNA was extracted according to the manufacturer's instructions of the miRNeasy mini kit (Qiagen). Equal amounts of isolated RNA (10ng) were subjected to one-step real-time PCR using the One-Step SYBR PrimeScript RT-PCR kit (Takeda Bio) on an HT7900 system in order to determine the percentage of RNA retrieval. DNA, on the other hand, was eluted twice using 150µL DNA elution buffer (50mM NaHCO3, 1% SDS, 200mM NaCl, 1mM PMSF, 100µg/mL RNase A (Thermo Fisher Scientific), 100U/mL RNase H (New England BioLabs)) and incubated for 30min at 37°C with shaking. DNA eluates were combined and incubated with 300µg proteinase K for 45min at 50°C. DNA was purified using phenol:chloroform:isoamyl (25:24:1, Thermo Fisher Scientific) followed by the MiniElute PCR purification kit (Qiagen). The quality of DNA was assessed using a 2100 Bioanalyzer instrument (Agilent). Finally, DNA was used to prepare libraries using the ThruPLEX DNA-Seq kit (Rubicon), and sequencing was performed according to the manufacturer's protocol on a HiSeq system (Illumina).

ChIRP-Seq analysis
Analysis including alignment and peak calling was performed following the published ChIRP-Seq pipeline from the Chang Lab at Stanford University 54 with minor modifications. Briefly, raw reads were firstly trimmed using SolexaQA (ver. 3.1.7.1) 100 with -dynamictrim and -lengthsort. Trimmed reads were then uniquely mapped to human reference genome hg38 using bowtie (ver. 1.1.1) 101 . Mapping parameters were -m 1 -chunkmbs 1024 -p 6. Peaks against input were called using MACS 2.0 (ver. (n = 124,047) and total genes from ChIRP and ASOKD data (n = 13,127) were respectively used. The circular plot of the identified gene targets was generated using Circos (ver. 0.69-7) 104 .
The empty pCDH plasmid was used as a negative control. PCR cloning primers are listed in Supplementary Table 11. Lentivirus production for each vector was carried out as described above.

Establishment of LECs overexpressing LESR2 and SEMA3C
1.2x10 5 LECs per well were seeded into a 6-well plate and cultured overnight. LECs were then infected with medium containing viruses overexpressing LESR2 and SEMA3C diluted at a 10 MOI and 10µg/mL polybrene (hexadimethrine bromide). Plates were then sealed with parafilm and centrifuged at 1,250rpm for 1.5h at room temperature. After 16-24h, the virus-containing medium was changed. 24h later, infected LECs were subjected to antibiotic selection using puromycin at a 1µg/mL concentration (Thermo Fisher Scientific). Once plates were confluent, infected LECs were split into 10cm dishes at 3x10 5 cells per dish and further selected with 2-5µg/mL puromycin until full confluence was reached. Finally, infected LECs were aliquoted and frozen down for further experiments. To check RNA levels, approx. 70'000 infected LECs were lysed, total RNA was extracted, and cDNA synthesis and qPCR were performed as described above. To check protein levels, a confluent 6cm dish of infected LECs was lysed in 350µL lysis buffer (25mM HEPES, 5mM EDTA, 1% triton-X, 150mM NaCl, 10% glycerol, complete protease inhibitor cocktail). The lysate was then centrifuged at 13,000rpm for 20min at 4°C, and the supernatant was collected in a 1.5mL Eppendorf tube. SEMA3C protein levels were checked by western blot as described above using rabbit anti-human SEMA3C antibody (1:1000, Thermo Fisher Scientific). Beta-actin (housekeeping gene) was detected with a rabbit anti-beta-actin antibody (1:5000, Abcam) (Supplementary Figure 6b).

Rescue of cell growth phenotype experiments
For LESR2, 5x10 5 LECs infected with pCDH-EV or pCDH-LESR2 were seeded into 10cm dishes and cultured overnight in starvation medium (EBM supplemented with 1% FBS). The next day, infected LECs were transfected with scramble ASO or LESR2-ASO2 for 24h as described above.
For KLF4, a consecutive transfection of ASOs and siRNAs was performed. siRNA sequences are listed in Supplementary Table 11. 350'000 LECs were seeded into 10cm dishes and cultured overnight in starvation medium. The next day, LECs were transfected with scramble control ASO or LESR2-ASO2, as described above. 24h post-transfection with ASOs, the medium was changed, and LECs were treated for additional 24h with 20nM of high GC scramble control siRNA or two siRNAs targeting KLF4 (Thermo Fisher Scientific) and 32µL Lipofectamine RNAiMAX previously mixed in 800µL Opti-MEM.
At the end of the experiment, LECs were detached and collected. 1.5x10 5 LECs per replicate were then transferred in a 96 U-bottom plate (Greiner bio-one). Aliquots of approx. 1x10 5 LECs were lysed and subjected to qPCR, as described above. Cell cycle analysis using flow cytometry was performed as described above. To detect Ki-67, an eFluor 450-conjugated rat anti-human Ki-67 antibody (clone: SolA15, ebioscience) was used at a 1:200 dilution.

Rescue of cell migration phenotype experiments
2.5x10 5 LECs infected with pCDH-empty vector (EV), pCDH-LESR2, and pCDH-SEMA3C were seeded into 6cm dishes and cultured for at least 6h. Infected LECs were then transfected with 20nM of scramble control ASO or LESR2-ASO2 and 8µL Lipofectamine RNAiMAX previously mixed in 800µL Opti-MEM according to the manufacturer's instructions. 24h post-transfection, LECs were detached, and 15,000 were seeded into a 96-well plate and cultured overnight. The next day, infected LECs were pre-treated with proliferation inhibitor, scratched, and monitored for cell migration as described above.

Biotin-RNA pull-down followed by mass spectrometry
Biotin-RNA pull-down experiments were performed as previously described 57  Lysed nuclei were transferred into a 7mL Dounce homogenizer (Kimble) and sheared mechanically using 30-40 strokes with pestle B. Next, the nuclear lysate was transferred to a fresh tube and sonicated 2x 30s at a high intensity (50% cycle and 90% power) using a Sonopuls HD2070 (Bandelin). 10U/mL DNase I (Thermo Fisher Scientific) were subsequently added to the nuclear lysate and incubated for 30min at 4°C while rotating. The nuclear lysate was further sonicated for another 2x 30s at high intensity. The nuclear lysate was centrifuged at 13,000rpm for 10min at 4°C.
Finally, the supernatant was collected into a fresh tube, and glycerol was added to reach a 10% final concentration. Resulting clear nuclear lysate was flash-frozen in liquid nitrogen and stored at -80°C for later use. Nuclear fractionation was checked after performing western blot analysis, as described above, of GAPDH (cytoplasmic protein) and Histone H3 (nuclear protein) using rabbit anti-GAPDH antibody (1:5000, Sigma) and rabbit anti-histone H3 antibody (1:10000, Sigma) antibodies (Supplementary Figure 7a).
To produce biotinylated RNA, full-length LESR2-1, determined by RACE (Supplementary Table 7), and antisense strand were cloned into pcDNA3.1 backbone. Both transcripts were amplified by PCR using Phusion high-fidelity DNA polymerase (New England BioLabs), according to the manufacturer's protocol. Amplified fragments were digested with BamHI and XbaI restriction enzymes (New England BioLabs) overnight at 37°C. After gel purification, digested fragments were cloned into a linearized pcDNA3.1 backbone (from 64599, Addgene) using T4 DNA ligase for 20min at room temperature.
Cloned vectors were transformed into one shot TOP10 chemically competent cells, according to the manufacturer's instructions. Plasmids were isolated using the Nucleospin Plasmid kit (Machery Nagel). Sequences of inserted fragments were checked by Sanger sequencing (Microsynth). PCR cloning primers are listed in Supplementary Table 11. Subsequently, both transcripts were biotinlabeled after in vitro transcription from 1µg linearized pcDNA3.1-LESR2-1 and pcDNA3.1-LESR2-1antisense plasmids for 1h at 37°C using Ampliscribe T7-flash biotin-RNA kit (Lucigen). Biotinylated LESR2 sense and antisense RNA were then treated with RNase-free DNase I for additional 15min at 37°C. Both biotinylated RNAs were purified by ammonium acetate precipitation, as described by the manufacturer. After determining the concentration using Nanodrop 1000, the integrities of sense and antisense LESR2 transcripts were tested by gel electrophoresis (Supplementary Figure 7b).
To perform RNA pull-down, 150µL Dynabeads M-270 streptavidin magnetic beads were washed twice with RNA pull-down buffer. For each condition, 60µL washed beads were then incubated with 1.5mg nuclear lysate for 30min at 4°C. During nuclear pre-clearing, 100pmol per condition of biotinylated RNAs were denatured by heating to 65°C for 10min and cooled down slowly to 4°C. Pre-cleared nuclear extract was further diluted to 2mL using RNA pull-down buffer and incubated with 100pmol biotinylated RNA for 1h at 4°C on a rotatory shaker. Next, 60µL washed streptavidin magnetic beads were added and further incubated for 45min at 4°C. Beads were carefully washed five times in RNA pull-down buffer. Bound proteins were finally eluted twice by adding 3mM biotin in PBS (Ambion) to the beads and incubating them for 20min at room temperature and for 10min at 65°C. Eluted proteins were subjected to protein identification by mass spectrometry at the Functional Genomics Center Zurich (FGCZ). Proteins were pelleted by TCA precipitation using 20% TCA. Protein pellets were washed twice with cold acetone. Dry pellets were then dissolved in 45µL trypsin buffer (10mM Tris, 2mM CaCl2, pH 8.2) plus 5µL trypsin (100ng/µL in 10mM HCl) and 1.5µL 1M Tris pH 8.2. After microwaving for 30min at 60°C, dried samples were dissolved in 0.1% formic acid. Digested peptides were finally analyzed by LC/MS.

Analysis of RNA-pull down data
Database searches were performed using Mascot software (Matrix Science). Search results were then loaded into Scaffold software (ver. 4.8.7) to perform statistical analysis. Only proteins with 1% protein FDR, a minimum of 2 peptides per protein, 0.1% peptide FDR, and present in both LESR2 replicates but not in the antisense control were considered. PPI network for the proteins identified by RNA-biotin pull-down was generated using the STRING web tool (https://string-db.org/cgi/input.pl) 58 .
The human PPI database was used for the analysis, while default values were used for the rest of the parameters. The identified proteins are listed in Supplementary Table 10.

Native RNA immunoprecipitation followed by qPCR
RNA immunoprecipitation experiments were performed as previously described 105 with minor modifications. To prepare nuclear lysate, 40 million LECs per replicate were collected and washed once with DPBS. The cell pellet was resuspended in 40mL consisting of 8mL DPBS, 24mL RNasefree H2O, and 8mL nuclear isolation buffer (1.28M Sucrose, 40mM Tris-HCl pH 7.5, 20mM MgCl2, 4% triton X-100). After mixing by inversion, LECs were let stand on ice for 20min with occasional mixing. Cells were then centrifuged at 1,668rpm for 15min at 4°C. Nuclear pellet was resuspended in 2mL of RIP buffer (150mM KCl, 25mM Tris-HCl pH 7.5, 5mM EDTA pH 8, 0.5% NP-40, 0.5mM DTT, complete protease inhibitor cocktail, 100U/mL Ribolock RNase inhibitor) and incubated for 5min on ice. Lysed nuclei were transferred into a 7mL Dounce homogenizer (Kimble) and sheared mechanically using 30-40 strokes with pestle B. Next, the nuclear lysate was transferred to a fresh tube and centrifuged at 13,000rpm for 10min at 4°C. Finally, the supernatant was collected into a fresh tube, and glycerol was added to reach a 10% final concentration. The resulting clear nuclear lysate was snap-frozen in liquid nitrogen and stored at -80°C for later use. For each replicate, 1mg nuclear lysate in 1.5mL RIP Buffer was incubated with 2.5µg rabbit anti-human RBBP7 antibody (Cell Signaling) or rabbit IgG isotype control antibody (I5006, Sigma) overnight at 4°C with gentle rotation.
Two aliquots of 15µL nuclear lysate were used as input RNA and protein. Next, 50µL washed Dynabeads protein A magnetic beads (Thermo Fisher Scientific) were added to each sample and incubated for 2h at 4°C. After five washes with 500µL RIP buffer, each sample's beads were resuspended in 500µL RIP buffer. 25µL beads of each sample were then transferred to a fresh tube and subjected to western blot analysis as described above using the same RBBP7 antibody at a dilution of 1:1000 (Figure 7c). In order to prevent masking from denatured IgG heavy chains, we detected the primary antibody with a mouse anti-rabbit IgG conformation-specific secondary antibody (Cell Signaling). The rest of the sample's beads were resuspended in 700µL Qiazol (Qiagen), and RNA was extracted using the miRNeasy mini kit, according to the manufacturer's instructions.
Isolated RNA was then subjected to cDNA synthesis and qPCR, as described above. LESR2 Ct values were normalized to the housekeeping gene GAPDH.

Statistical analysis
All statistical analyses were performed using GraphPad Prism software (ver. 7.0.0). P-values were calculated after performing ordinary and RM two-or one-way ANOVA with Dunnett's correction, paired or unpaired Student's t-test as indicated. Statistical significance was determined when P < 0.05. If not alternatively specified, all error bars represent mean values with SD.

Data availability statement
Further information and requests for resources and reagents should be directed to and will be fulfilled by Michael Detmar (michael.detmar@pharma.ethz.ch). All software used in this study are specified in the related method paragraph with the corresponding version. All unprocessed sequencing data are going to be deposited in DDBJ DRA (https://www.ddbj.nig.ac.jp/dra/index-e.html) and on the FANTOM consortium website (https://fantom.gsc.riken.jp/6/). LC/MS data are available at the Functional Genomic Center Zurich (Yolanda Auchli, yolanda.auchli@fgcz.ethz.ch; ID:11594).

Acknowledgments
This study was financially supported by the ETH Zurich (grant ETH-24 171), the Swiss National Science Foundation (grant 310030_166490), and the European Research Council (advance grant LYVICAM). We acknowledge Jeannette Scholl for her support in performing the smRNAFISH expression analysis, as well as Dr. Raffaella Santoro and Dominik Bär, University of Zurich, for crucial advice on the RNA pull-down assay. Finally, we also thank all the members of the FANTOM6 project for fruitful discussions and support throughout the project.

Declaration of Interests
The authors declare no competing interests.  Table   1.      Supplementary   Table 5.     Data are displayed as mean + SD (n = 10 in a, f, and k; n = 5 in b; n = 3 in h, i, and j; n = 2 in d).

Figure Legends
Percentages represent LESR2 knockdown efficiencies after the experiments. **P < 0.01, ****P < 0.0001 using ordinary one-way (for a, b, d, and j) and two-way (for a, b, f, and k) ANOVA with Dunnett's multiple comparisons test against control ASO/sgRNA or LESR2-ASO2 -control siRNA. All displayed in vitro assays were performed in neonatal LECs derived from the same donor.       Data are displayed as mean + SD (n = 3 in b-d; n = 10 in e). Percentages represent the knockdown efficiencies of LESR2 after the experiments. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, ns not significant using RM one-way (for b, d), ordinary one-way (for c), and ordinary two-way (for e) ANOVA with Dunnett's multiple comparisons test against control ASO -control siRNA (b), LESR2-ASO2 -control siRNA (c), pCDH-EV -control ASO (d), and pCDH-EV-ASO2 (e). All displayed in vitro assays were performed in neonatal LECs derived from the same donor. (c) Representative western blot after RNA immunoprecipitation for RBBP7 followed by qPCR for LESR2 in neonatal LECs. To prevent masking from IgG heavy chain, a conformation-specific IgG secondary antibody was used. Uncropped western blot image is shown in Supplementary Figure 8. (d) LESR2 enrichment displayed as FC against IgG Control after RNA immunoprecipitation for RBBP7 in neonatal LECs derived from the same donor. LESR2 qPCR levels were normalized to the housekeeping gene GAPDH. Bars represent mean ± SD (n = 3). *P < 0.05 using unpaired Student's ttest against IgG Control.           Cell growth profile and cell growth rate (48h) after LESR2-ASOKD