Key molecular alterations in endothelial cells in human glioblastoma uncovered through single-cell RNA sequencing

Passage of systemically delivered pharmacological agents into the brain is largely blocked by the blood-brain-barrier (BBB), an organotypic specialization of brain endothelial cells (ECs). Tumor vessels in glioblastoma (GBM), the most common malignant brain tumor in humans, are abnormally permeable, but this phenotype is heterogeneous and may differ between the tumor’s center and invasive front. Here, through single-cell RNA sequencing (scRNA-seq) of freshly isolated ECs from human glioblastoma and paired tumor peripheral tissues, we have constructed a molecular atlas of human brain ECs providing unprecedented molecular insight into the heterogeneity of the human BBB and its molecular alteration in glioblastoma. We identified 5 distinct EC phenotypes representing different states of EC activation and BBB impairment, and associated with different anatomical locations within and around the tumor. This unique data resource provides key information for designing rational therapeutic regimens and optimizing drug delivery.


Introduction
Glioblastoma (GBM) is the most aggressive and lethal type of brain tumor, and the overall survival of the patients has not been improved substantially over the past 30 years (1). The poor prognosis for GBM patients is at least partially attributed to the extremely limited therapeutic options that are available after surgery and radiotherapy.
Despite tremendous basic and clinical research in the GBM field, only 4 drugs have been approved by the FDA for GBM treatment: temozolomide, lomustine, carmustine and bevacizumab (https://www.cancer.gov/about-cancer/treatment/drugs/brain) (2).
However, the effectiveness of these and other anti-tumor compounds, including the vast majority of low-molecular-weight chemotherapeutic drugs are thwarted by the unique features of the brain's blood vessels, known as blood-brain barrier (BBB), which tightly regulate the homeostasis of CNS (3,4). Failure to cross the brain vessel wall is likely a major contributor to the negative outcome of clinical trials for many blood-borne drugs (2).
Endothelial cells (ECs) are the key cellular component of the BBB. Brain ECs establish continuous complexes of tight and adherens junctions along EC-EC contacts providing a tight and size-selective barrier. They further express distinct sets of influx and efflux transporters, several of the latter with the ability to bind and efflux a broad variety of small lipophilic xenobiotic compounds and drugs (3). Brain EC also display very low levels of vesicular transcytosis further limiting passage of blood-born water-soluble molecules of all sizes (3). Molecular alteration of ECs has been observed in patients and animal models of many brain diseases, including stroke, multiple sclerosis, traumatic brain injury and GBM (5,6). By analysis of bulk mRNA isolated from brain ECs, we and others have shown that the abnormal vessels in GBM are associated with a distinct gene signature (7)(8)(9). However, the more precise characteristics of endothelial gene expression in GBM, including possible heterogeneity that cannot be resolved using bulk methods, is still poorly understood. Vascular targeting has been tried in order to prune and normalize tumor vessels and thereby improve drug delivery in GBM (10)(11)(12)(13). However, treatment of patient with bevacizumab to block VEGF signaling did not improve overall survival in unselected GBM patients, and also led to a decrease in temozolomide delivery (14).
To understand whether and to which extent ECs in GBM are phenotypically different from control brain ECs at single cell level, and to identify novel angiogenic pathways in GBM, we have performed scRNA-seq on ECs from human GBM tumor core and tumor peripheral tissue to characterize the heterogeneity of gene expression signatures of different ECs. We have identified five EC clusters and found that they were associated with distinct anatomical localizations and molecular phenotypes. Notably, the expression of many BBB-specific transporter mRNAs in GBM ECs was found to be retained to a larger extent than anticipated, and heterogeneous rather than uniformly downregulated or absent, together indicating that the BBB phenotype is significantly preserved even among the most abnormal tumor ECs in GBM. To our knowledge, this is the first scRNAseq-based molecular atlas of the human BBB and its perturbation in GBM. This unique resource will provide deeper insights into understanding the characteristic of endothelial cells in CNS malignancies, and information for designing rational therapeutic regimens and optimizing drug delivery.

scRNA-seq and cell type identification of glioblastoma and peripheral tissue
We isolated cells from two distinct regions in each of the surgically removed GBM tissue from 4 GBM patients: 1) the tumor core and 2) immediately neighboring brain tissue ( Figure 1A, Table S1). All the patients were diagnosed as new primary GBM and had not received chemo-or radio-therapy before surgery. Histopathological analysis confirmed a typical GBM picture, including a dense cellularity and abnormal blood vessel profiles within the tumor core, whereas the surrounding tissue showed an overall normal brain histological picture ( Figure 1B, 1C). GBMs are invariably recurrent in the perioperative region, and we therefore assume that the surrounding tissues contained tumor cells, albeit at low abundance and without obvious impact in the tissue histology.
Importantly, this normality included the vascular morphology, which was reminiscent of that of the normal brain ( Figure 1B, 1C). Magnetic-activated cell sorting (MACS) was used to generate CD31 + cell suspensions, which were subsequently subject to scRNA-seq using a 10X genomics protocol ( Figure 1A). CD31 is an endothelial antigen, which in human is expressed abundantly also on hematopoietic cells. Due to the limited amount of tissue available (average 1.5 g/sample), additional CD45 negative selection with the aim of removing leukocytes and enriching for endothelial cells failed to yield sufficient numbers of cells for an optimal 10x experiment. Thus, CD45 negative selection was not performed. After quality filtering, we captured 97,584 cells containing on average transcript reads from 2,224 genes per cell (Table S2).
The cells were clustered using the Seurat package (15) and annotated according to expression of canonical cell class markers ( Figure 1D, 1E, Table S3). In addition to ECs  Figure 1D, 1E, S1A-S1D).

Single-cell atlas of EC phenotypes in human brain and glioblastoma
The ECs were in silico selected based on canonical marker expression (CLDN5, VWF, CD34) and further analyzed. After removal of heterotypic cell doublets, re-clustering these revealed 5 EC clusters (cluster 1-5) (Figure 2A). Two of the clusters (cluster 1 and cluster 4) were almost entirely originating from brain tissue with non-malignant morphology surrounding the GBM tumor mass, whereas the remaining three clusters (cluster 2, cluster 3 and cluster 5) were largely derived from the tumor core ( Table S5). Interestingly, Co2 ECs also expressed high level of genes identified as tip cell marker in previous study, including FABP1A, CALM1, MARCHSL1 and GNG11 (16). We also identified 2 EC clusters that exhibit immune- This cluster was mainly derived from brain tissue surrounding the tumor core ( Figure   3D, Table S5). Cluster 5 (tumor core endothelial cell type III: Co3), were mainly derived from tumor core and were characterized by upregulation of immune-activated genes including IL1B, ACKR1, SELE and VACM1, which are associated with inflammation and immune cell recruitment ( Figure 3E, Table S5). Both these immuneactivated EC types (Pe2 and Co3) predominantly consisted of cells originating from individual patients ( Figure 2D), suggesting that they may represent individual differences in genetic background or systemic inflammation.

Tumor ECs acquire similar phenotypes across different tumor types
To explore whether tumor ECs have similar phenotype across different tumor types, we analyzed our data set together with a recent publicly available scRNA-seq dataset of ECs from human lung cancer and non-malignant lung tissue using Jaccard similarity analysis ( Figure 3F) (16). This analysis scored similarity of the top marker genes for all EC subclusters. As expected, Pe1 ECs from peripheral tissue exhibited a specialized phenotype, and were clearly distinct from other EC sub-populations ( Figure 3F).
Interestingly, we found that Pe2 ECs in peripheral tissue express similar markers and gene sets, including antigen presentation and immune/inflammatory response, to socalled "scavenging capillary ECs" from lung ( Figure 3E, S3A, S3B, Table S6). Notably, in tumors, Co1 and Co2 ECs in GBM closely resemble the phenotypes of tip cells and activated PCV in lung cancer, respectively ( Figure 3F, S3C-S3F). Several congruent genes were identified as top-50 markers for both Co1 EC in GBM and tip EC in lung cancer ( Figure S3C, Table S6). Gene sets involved in extracellular matrix organization, angiogenesis and cell migration were commonly up-regulated in both EC subtypes ( Figure S3D, Table S6). Co2 EC in GBM and activated PCV in lung cancer shared common markers and gene sets associated with mRNA translation and ribosome assembly ( Figure S3E, S3F, Table S6). Taken together, these results suggest that although GBM and lung cancer are distinct tumors, they comprise endothelial cells exhibiting partially overlapping gene signatures.

Distinct anatomical localization of individual EC cluster elucidated using reference atlases
We next focused our attention on the phenotypes of Pe1 ECs, Co1 ECs and Co2 ECs, which were represented in multiple patients. We analyzed the association of the EC clusters with anatomical location using the Ivy GAP database, which contains transcription profiles of human GBM anatomic regions including leading edge, infiltrating tumor region, cellular tumor region, microvascular proliferation and pseudopalisading necrosis region (http://glioblastoma.alleninstitute.org/) (17). We found that normalized top-50 markers of Co1 ECs were enriched in microvascular proliferation regions and down-regulated in leading edge region ( Figure 3G). In contrast, Pe1 EC markers, including SLCO1A2, ANXA3, TSC22D1, were enriched in leading edge and infiltrating tumor regions ( Figure 3G). These results confirm the distinct anatomical localization of ECs with different phenotypes using an independent dataset.

ECs in GBM are distinct from ECs in peripheral brain tissue
In order to explore which transcription factors may regulate the distinct phenotype of ECs in tumor periphery and tumor core, we employed Single-Cell Regulatory Network Inference and Clustering (SCENIC) (18) to evaluate the activated transcription factors in Pe1 and tumor core ECs including both Co1 and Co2 (Table S7). This identified several activated transcription factors in ECs in tumor periphery and tumor core respectively ( Figure 4A, Table S7). Notably, SOX4 and ETS1, which were identified as activated transcription factors in ECs in the tumor core, were also up-regulated in ECs in the tumor core ( Figure 4C). Interestingly, elevated ETS1 expression was also observed in vasculature of IDH-wild type lower grade gliomas (LGGs) compared to vasculature of IDH-mutated LGGs in our previous study (9). SOX4 has been suggested to promote tumor angiogenesis through CXCR4 and endothelin-1 in breast and hepatocellular tumors respectively in two recent studies (19,20).
To evaluate global metabolic alteration in ECs, we performed GSVA with metabolic gene sets (21) for ECs in tumor peripheral tissue and tumor core. This analysis revealed that ECs in tumor core displayed an upregulation of glycolysis, citric acid cycle and oxidative phosphorylation gene expression signatures ( Figure S4A), similar to what was previously observed in ECs from human lung cancer as compared to non-malignant lung (22). Considering that ECs are highly glycolysis-addicted cells, high glycolysis in tumor ECs may reflect high demand of energy requirements for angiogenesis in the tumor microenvironment. Indeed, therapeutically inhibiting glycolysis could inhibit tumor angiogenesis and normalize tumor vessels (23). Up-regulation of oxidative phosphorylation and the tricarboxylic acid (TCA) cycle are noteworthy, and this observation is in line with a recent study showing that the oxidative phosphorylation is necessary for angiogenesis (24). Inhibition of oxidative phosphorylation by ablation of a subunit in respiratory chain complex III in ECs leads to diminished EC proliferation and impairment in retinal and tumor angiogenesis, accompanied by decreased amino acid levels but without affecting on anabolism or nucleotide levels (24). To further explore which metabolic genes and gene sets were regulated in ECs in tumor, we performed differential gene expression analysis and constructed a map of pathways involved in central carbon metabolism ( Figure S4B). The results confirmed upregulation of genes in metabolic pathways supporting biomass synthesis including glycolysis, citrate cycle, oxidative phosphorylation and nucleotide synthesis and downregulation of genes in glutamate metabolism in ECs in tumor core ( Figure S4B).
Taken together, our result indicates that ECs in tumors are associated with an altered metabolic transcriptome signature characterized by increased expression of genes involved in glycolysis, citrate cycle and oxidative phosphorylation.

Identification of novel GBM endothelial markers and angiogenic regulators
In order to identify novel angiogenic regulators and to discover markers discriminating

ECs in GBM are associated with partially intact BBB phenotype characterized by downregulation of transporter genes and upregulation of transcytosis gene
Emerging studies suggest that the BBB is disrupted during tumor progression (reviewed in (6)). However, whether and to what extent BBB alterations in GBM ECs can be documented at the single-cell level remains unknown. We therefore compared our data to a recently reported BBB dysfunction module, which compiles a set of genes that are upregulated in ECs in at least 3 out of 4 mouse brain disease models where the BBB has been disrupted, including stroke, multiple sclerosis, traumatic brain injury and seizure (5,25). BBB core module comprises a group of genes enriched in brain ECs KLF2, a key transcription factor orchestrating a network of genes that promotes EC quiescence in response to flow (29), is one of the top-10 enriched genes in brain EC cluster. GLUT1, encoded by SLC2A1, is highly expressed in BBB ECs, and facilitates glucose transport over BBB (25). Depletion of GLUT1 in adult brain ECs leads to activation of inflammatory and extracellular matrix-related gene sets (25).
We also identified ECs with angiogenic phenotype (cluster 2: Co1) that expressed high It is common belief that the BBB is disrupted during tumor progression in the brain.
Our scRNA-seq data indeed show that tumor ECs in GBM have a partially intact BBB phenotype, characterized by downregulation of transporter genes, whereas the expression of junctional molecules remained normal or was increased. ECs in GBM express high level of plasma lemma vesicle-associated protein (PLVAP), which is a vascular marker of BBB disruption. PLVAP expression is absent in brain vasculature, with the exception of the choroid plexus and circumventricular organs where the endothelium is fenestrated to allow filtration of plasma for CSF production and passage of hormones (31). However, PLVAP is induced in pathological conditions in the brain and associated with vascular leakage (28). PLVAP is a key regulator of vascular permeability and promotes transcytosis in ECs by forming the diaphragms of caveolae, fenestrae and trans-endothelial channels (28). Vascular leakage is a hallmark of GBMs and may be induced by two main pathways: (1) increased paracellular transport by altering the tight junctions between ECs and/or (2) increased transcellular transport by altering vesicular transcytosis. In contrast to upregulation of PLVAP in Co1 ECs, the expression of tight junction molecules including CLDN5 was similar between ECs in periphery and tumor core, highlighting an important role of transcellular pathways for BBB breakdown and the edema formation observed in GBM.
Several key BBB transporters, including SLC2A1, ABCG2, ABCB1, SLCO1A2 and ATP10A, were highly expressed in the ECs of the brain tissue surrounding the GBM core. P-glycoprotein (P-gp; encoded by ABCB1) and breast cancer resistance protein This resource may provide vital information of relevance for drug delivery and intratumoral distribution in GBM, and facilitate the design of rational therapeutic regimens.

Isolation of CD31 + endothelial cells from human glioblastoma and peritumoral tissue.
We collected surgical tissues from 4 GBM patients to isolate ECs for scRNA-seq. The GBM tumors were located in the right parietal lobe for patient 1; the left temporal lobe for patient 2; the right occipital lobe for patient 3 and the right frontal lobe for patient 4 (Table S1). From each patient, we collected two separate tissue samples: one originating from the tumor core and another from the peri-tumoral space. For each sample, tumor core and peritumor tissue were processed separately. Tissue samples were transported to the research facility immediately from the operating room in order to start sample dissociation within 2 hours of resection. Tissue samples were mechanically dissociated and then processed into single cell suspension using tumor dissociation kit (130-095-929, Miltenyi Biotec) for tumor core tissue; using adult brain dissociation kit (130-107-677, Miltenyi Biotec) for peritumor tissue. Single cell suspensions were run through debris removal solutions to remove the myelin debris and then proceeded to red blood cell removal according to manufacturer's specifications. CD31 + selection for endothelial cell enrichment was performed using Dynabeads (11155D, Invitrogen). For further processing for scRNA-seq, the samples were resuspended in DPBS containing 0.04% bovine serum albumin. The number of cells and fractions of live cells in suspension was counted and the volume of suspension containing the required number of live cells were used for scRNA-seq as described below.
Single-cell Droplet-Based single cell sequencing, quality control (QC) and data processing scRNA-seq libraries were prepared using Chromium Single Cell Reagent Kit (10X Genomics). The libraries were then pooled and sequenced on NovaSeq 6000 (Illumina), using NovaSeq Control Software v1.6.0. The raw sequence data were processed using Cell Ranger software (v.3.0.1). The reads were aligned to human genome GRCh38 and gene counts matrix was generated for each sample. The raw counts data were then loaded into Seurat package (v3.1.1) for quality control, filtering, normalization, UMAP visualization and clustering (15). The cells which have mitochondrial genes greater than 10% or have less than 200 detected genes were filtered out. A scale factor of 10,000 was used to normalize all the remaining cells. To correct for the batch effect between different samples, the rpca method in the Seurat package was applied to integrate all the dataset. The genes enriched in each cluster were identified using FindAllMarkers function. It applies a Wilcoxon Rank Sum test and then performs multiple test correction using Bonferroni method. The multiple test corrected p value < 0.05 was used as cutoff for significance.

Similarity analysis for different tumor ECs
To illustrate the similarity of the different tumor ECs groups from our study with previous published lung tumor EC subtypes (16), we applied Jaccard similarity analysis.
The top-50 enriched genes from each cluster were compared and a pair-wise Jaccard similarity coefficients matrix was calculated. The result matrix was then visualized in 2-D using the Classical Multidimensional Scaling method in R (version 3.6.1)

Ivy GAP database
The transcriptome data from different human GBM anatomic regions, including leading edge, infiltrating tumor region, cellular tumor region, microvascular proliferation and pseudopalisading necrosis region, was obtained from Ivy GAP database (http://glioblastoma.alleninstitute.org). Due to heterogeneity of vascular abundance in different anatomical localization, a direct comparison of different EC cluster markers is inappropriate. Therefore, the original marker gene expression values were normalized by a microvascular score, which was calculated by vascular enriched genes to estimate the abundance of vasculature, as described previously (9). The normalized expression of the marker genes for different EC clusters in distinct anatomical localization were shown by heatmap as Figure 3G.

Gene set variation analysis (GSVA)
To identify the gene sets with significant changes in the tumor clusters, we applied GSVA analysis using GSVA package (version 1.32.0). The metabolic gene sets obtained from a published study (21), were tested using the R limma package (version 3.40.6).
Multiple-test was corrected using the Benjamini-Hochberg method and the gene sets with corrected p value < 0.05 were identified as significant.

SCENIC analysis
To identify the transcription factors that regulate the tumor EC clusters, SCENIC analysis was performed using the RcisTarget package (version 1.4.0) using default settings. It identifies overrepresented transcription factor binding motifs on a gene list, which were then annotated to transcription factors.

Core BBB and BBB dysfunction module
Core BBB module gene set was generated according to the previous study which compared transcriptome of mouse brain ECs to ECs in peripheral organs including kidney, lung, heart and liver (5). Core BBB module gene set comprise 162 genes which were selected based on following criteria: 1. at least 100 CPM (counts per million mapped reads) detected in brain ECs; 2. more than 2-fold (and P＜0.05) in brain ECs compared to ECs from other 4 organs individually; 3. with expression in brain ECs no less than brain vasculature (exclude contamination from mural cells). 155 out of the 162 genes have human homologs, and 147 of them were detected in our data set.
BBB dysfunction module gene set comprise 136 genes which are upregulated in ECs in at least 3 out of 4 mouse diseased models including stroke, multiple sclerosis, traumatic brain injury and seizure in previous study (5). 131 genes have human homologs, and 128 of them were found in our data set.

Identification markers for control and GBM vasculature
Direct comparison of ECs to other CD31 positive cell types yield to 374 EC enriched genes (Bonferroni-corrected p value < 0.05). The expression of 374 EC enriched genes were compared between ECs in periphery and tumor core subsequently, listed in Table   S8.

HE staining analysis
Patient samples were fixed in 4% formalin, embedded in paraffin, followed by section and staining with hematoxylin and eosin (H&E). Then the sections were examined for the presence of tumors.

Immunohistochemical staining and quantification
Immunohistochemistry was performed on 6-μm sections of formalin-fixed, paraffin- (NBP2-44448, Novus) overnight at 4℃, followed by incubation with the secondary antibody and nuclear staining with Hoechst 33258 (Sigma-Aldrich). The slides were then mounted with Fluoromount (Sigma-Aldrich). All the images were acquired by Axio Imager upright microscope (Zeiss, Germany).

Data availability
The sc-RNAseq raw sequencing data and also processed counts data are available in the NCBI Gene Expression Omnibus under accession number GSE162631.

Statistics
The immunohistochemistry stainings were quantified and analyzed using the two-tailed Wilcoxon test. A P value of less than 0.05 was considered significant.

Study approval
Samples were obtained following informed consent, under the auspices of the Tangdu Hospital of Fourth Military Medical University (THFMMU, 2019-0166). Ethical permit of using patient samples was granted by ethics committee of Shaanxi Normal University, and informed consents were obtained from all subjects (Table S1). Any information that might disclose the identity of the subjects has been omitted.