Single-cell transcriptomic profiling provides insights into retinal endothelial barrier properties

Isolation and culture of retinal microvascular endothelial cells

RMECs were isolated from bovine retinas using protocols well established in our laboratory that have been demonstrated previously to consistently generate highly pure endothelial cultures with minimal contamination [32]. Bovine eyes were transported from a local abattoir on ice. Ten eyes were used per RMEC isolation. The retinas were removed, washed free of RPE, and homogenized. Microvessel fragments were trapped on an 85 µm filter. The microvessels were subsequently digested using an enzyme cocktail of Pronase E from Streptomyces griseus (Cat No 1074330001, Sigma-Aldrich, Gillingham, UK), DNase1 (Cat No LS002139, Worthington Biochemicals, Lakewood, NJ) and Collagenase from Clostridium histolyticum (Cat No C9891, Sigma-Aldrich, Gillingham, UK) and passed through a 53 μm filter. They were then centrifuged (800 ×g for 10 min) and the pellet resuspended in Dulbecco’s modified Eagle’s medium (DMEM) containing 0.1 mg/ml Primocin (InvivoGen; San Diego, CA), 1.4 ng/ml insulin (Sigma-Aldrich), 5 µg/ml heparin (Sigma-Aldrich), and 20% (v/v) porcine serum (Sigma-Aldrich). This mixture was seeded into a 25 cm², 1% gelatin-coated culture flask (TPP, Trasadingen, Switzerland) and maintained at 37 °C in 5% CO₂. The culture medium was changed every 48 h for approximately 5–6 days before the cells were passaged. The RMECs were passaged in a 1:3 ratio at approximately 80% confluency by washing the cells with phosphate buffered saline NaCl (PBSA; 1X; 171 mM NaCl, 3.4 mM KCl, 10.1 mM Na₂HPO₄, 1.8 mM KH₂PO₄, pH 7.4) before TrypLE Express (Gibco, ThermoFisher Scientific, Leicestershire, UK) was added, and the flask was incubated for approximately 2 min to allow the cells to detach. The reaction was stopped with the addition of culture medium. From passage 1 onward, the cells were maintained in DMEM supplemented with 0.1 mg/ml Primocin, 1.4 ng/ml insulin, 5 µg/ml heparin, and 10% (v/v) porcine serum and used for experiments at passage 3.

Immunohistochemistry

Animal use conformed to the standards of the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research and was authorized by Queen’s University of Belfast Animal Welfare and Ethical Review Body (AWERB). Adult C57BL/6J mice were euthanized with CO₂ asphyxiation and cervical dislocation. The eyes were immediately enucleated, punctured, and fixed in 4% paraformaldehyde (PFA) in PBS (1X; 136.9 mM NaCl, 2.7 mM KCl, 8.1 mM Na₂HPO₄, 1.5 mM KH₂PO₄) for 1–2 h at room temperature. The retinas were dissected out and placed in permeabilization buffer (0.5% Triton X-100 with 5% normal donkey serum in PBS) overnight at 4 °C. Primary antibodies listed in Appendix 1 were selected based on their specificity toward the associated protein and prepared in permeabilization buffer with isolectin B4-biotin (1:200, Sigma-Aldrich). The retinas were incubated with primary antibodies and isolectin B4-biotin at 4 °C for 4 days. Following extensive washing in PBS for 8–9 h, the retinas were incubated with donkey anti-rabbit Alexa 488 and streptavidin Alexa 568 (both 1:200, ThermoFisher, Gillingham, UK) at 4 °C overnight. Nuclei were labeled with the far-red nuclear stain TOPRO3 (1:1,000, ThermoFisher) in PBS. The retinas were rinsed once with PBS and flatmounted on glass sides in Vectashield antifade solution (Vector Laboratories, Peterborough, UK). Images were acquired using a Leica SP5 confocal laser scanning microscope (Leica Geosystems; Heerbrugg, Switzerland; HCX PL APO x63 oil immersion lens). Each fluorescent channel was excited and captured sequentially to minimize bleed-through. Following acquisition, images were digitally segmented using the isolectin-B4 channel to isolate retinal capillaries from the surrounding retinal tissue. All secondary-only controls were negative for staining.

RNA isolation and RT–PCR

RNA was isolated from RMECs at passage 3 using the RNeasy Mini Kit (Qiagen, Manchester, UK) according to the manufacturer’s instructions. cDNA synthesis was performed using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, ThermoFisher Scientific). PCR was performed with REDExtract-N-Amp PCR Reaction mix (buffer, salts, dNTPs, Taq polymerase, REDTaq dye, and JumpStart Taq antibody; Merck, Darmstadt, Germany) in conjunction with gene-specific forward and reverse PCR primers (Appendix 2). Agarose gel electrophoresis was conducted at 200 V for 45 min on a 2.5% agarose gel made with 2.5 g agarose powder (Apollo Scientific, Manchester, UK) and 100 ml 1X lithium borate buffer (10 mM lithium acetate, 10 mM boric acid) [33,34]. The GeneRuler low-range DNA ladder (ThermoFisher Scientific) was used to identify band sizes.

Single-cell RNA library preparation

RMECs were grown to confluency at passage 3 and dissociated using TrypLE Express (Gibco, ThermoFisher Scientific). Cell concentration was determined using an Eve Automated Cell Counter (10 µl cells/10 µl Trypan blue, Nanoentek USA Inc., Loughborough, UK), and the suspension was diluted to 2 × 10⁵ cells/ml so that about 500 cells entered the cell inlet of a C1 Single-Cell Auto Prep integrated fluidic circuit (IFC) for mRNA Seq, 10–17 μm (Fluidigm, San Francisco, CA, Part No. 100–5760). Cell buoyancy was assessed according to the Fluidigm Single-Cell Preparation Guide (PN 100–7697). Cell capture sites were visualized using the Auto Find Fluidigm protocol on an EVOS FL Auto Imaging System (Thermo Fisher). cDNA synthesis and amplification were performed with a Clontech SMART-Seq v4 Ultra Low Input RNA Kit for Fluidigm C1 IFCs. cDNA harvested from the Fluidigm C1 was quantified using the AccuBlue High Sensitivity dsDNA Quantification Kit (Thermo Scientific) in combination with an Echo 525 Acoustic Dispenser (Labcyte, San Jose, CA). The Clontech modified Illumina® Nextera® XT DNA sample preparation protocol was also miniaturized using an Echo 525 Acoustic Dispenser to prepare single-cell mRNA sequencing libraries from cDNA acquired using the C1 system. cDNA was fragmented and tagged with adaptor sequences using a Nextera XT DNA Library Preparation Kit (Illumina). For each sample well (96 samples in a 384 Echo plate), 500 nl of Tagment DNA Buffer, 250 nl of Amplicon Tagment Mix, and 250 nl of sample (at 50 pg) were added to give a reaction volume of 1,000 nl. The plate was then sealed and centrifuged at 1,500 ×g for 30 s before incubation at 55 °C for 5 min and held at 10 °C. Once the sample reached 10 °C, 250 nl of Neutralize Tagment Buffer was added. Tagmented cDNA then underwent limited-cycle PCR amplification to add index adapters and sequences required for cluster formation. To each sample, 1.875 µl of Nextera PCR Master Mix, 1.875 µl of water, and 0.625 µl of i5 (S5XX) and i7 (N7XX) adapters were added to give a total reaction volume of 6 µl. The plate was sealed and centrifuged at 1,500 ×g for 1 min before incubation on a preheated thermocycler; 72 °C 3 min, 95 °C 30 s and then 12 cycles; 98 °C 10 s, 55 °C 30 s, 72 °C 30 min followed by 72 °C 5 min, and held at 10 °C. The libraries were purified using Agencourt AMPure XP Bead Clean-up (1.8X). The libraries were then pooled and normalized, and the molarity was quantified using the High Sensitivity NGS Fragment Analysis Kit on a Fragment Analyzer (Agilent Technologies).

Sequencing

Sequencing was performed by the Genomics Core Technology Unit (GCTU) unit at Queen’s University Belfast (Genomics) using the NextSeq® 500 system (Illumina).

Bioinformatic analyses

Fastq files were generated using bcl2fastq2 in the BaseSpace cloud platform (Illumina) and sequencing quality assessed using fastqc. Samples with no cell capture, >1 cell per capture site, or <500 pg/µl cDNA library were discarded before alignment. Reads were aligned to Bos taurus UMD3.1 genome using the Spliced Transcripts Alignment to a Reference (STAR) read mapper [35], and cells with <100,000 aligned reads were discarded. Alternative splicing within single cells was quantified using Outrigger [31], which uses the junction reads from the STAR output to quantify percent spliced in values (PSI /ψ) for each event within each single cell. The Integrated Genomics Viewer (IGV) [36] was used to visualize splicing and create sashimi plots of alternative splicing events. Gene counts were calculated using HTseq [37], and an expression matrix was created. Quality control (QC) was visualized at each stage using multiQC [38]. Seurat [39] is an R package that was used for QC, preprocessing, dimension reduction, clustering, and differential expression of scRNA-seq data. Seurat takes the raw expression matrix as input and was able to mitigate the effects of cell cycle heterogeneity in scRNA-seq data by calculating cell cycle phase scores and regressing them out of the data during preprocessing. The R-based Single-cell Analysis Toolkit for expression in R (Scater) [40] was used to generate violin plots of gene expression. The single-cell latent variable model (scLVM) [41] within the web-based single-cell transcriptomic analytical tool Automated Single Cell Analysis Pipeline (ASAP) [42] was used for normalization of gene expression, accounting for differences in sequencing depth and technical bias. Single-cell consensus clustering (SC3) [43] was used for unsupervised clustering of the filtered data set within ASAP [42]. The ideal number of clusters was determined by maximum silhouette value [43,44], which indicated the stability of each cluster. ASAP was also used for exploration and visualization of the data, specifically expression of genes related to the microvascular barrier.

The BAM files from STAR were analyzed using Monovar [45], a single nucleotide variant (SNV) detection and genotyping algorithm for single-cell sequencing data, to generate variant call files (VCFs). They were processed using VCFtools [46] in a matrix with individual cells in columns and SNVs in rows (labeled 0, 1, or 2 for variant absent, present on one allele, or present on both alleles, respectively). This matrix was used to cluster cells based on their specific SNV fingerprint into groups of cells originating from the same animal. The Cytoscape open source software platform [47] was used in conjunction with the GENEmania [48] and String [49] applications for network analyses and generation of the blood–retina barrier model. Protein functions and interactions were analyzed using the STRING database [50]. CorrelationCalculator was used to calculate and generate heatmaps of pairwise correlations between genes based on their expression pattern among individual cells [51]. Correlation results were imported into Cytoscape and MetScape used to view correlation networks [51,52].

Population structure of cultured RMECs

To assess the existence of distinct subgroups within RMECs, unsupervised clustering was performed to group the cells based on the similarity of their transcriptomes (i.e., all genes expressed within individual cells). SC3 [23] demonstrated that heterogeneity between RMECs was due principally to cell cycle variation, with proliferating cells forming a small distinct cluster (Appendix 4). All genes differentially expressed between the non-dividing and cycling cells are detailed in Appendix 5.

Use of the Seurat package [39] to regress out cell cycle effects confirmed that the RMECs otherwise formed a relatively homogeneous population (Figure 1). The RMEC cultures were prepared by pooling retinas from five cows; therefore, the genetic origin of each cell could be traced using its characteristic SNV profile. Cells from one individual animal did not form subclusters with similar gene expression (Figure 1, Appendix 4). This result suggests that the genetic background or individual differences in the donor cows’ phenotypes did not have a marked effect on an individual cell’s gene expression, but more cells would be required to draw any firm conclusions.

Arteriovenous marker genes have been well defined at the single-cell level for brain vasculature [53]. We observed predominant expression of genes highly expressed in brain capillary endothelial cells and limited expression of most arterial or venous markers, although wider expression of Vegfc and Nr2f2 (Figure 2).

Development of a BRB model: Structural junctional proteins

Maintenance of the iBRB is a key function of the retinal microvascular endothelium, and we used the genes expressed in RMECs as a predictor of the main proteins involved in this function. First, a subset of barrier genes listed by the Gene Ontology consortium (GO:0061028 - Establishment of endothelial barrier, GO:0005911 - cell–cell junction and GO:0071603 - endothelial cell–cell junction) that were expressed within the RMECs was created. Known interactions between the proteins encoded by these genes were added to create the model depicted in Figure 3 (the list of genes in the network is provided in Appendix 6, and the Cytoscape session file is available upon request). Several of the proteins in this model have not been previously shown in RMECs or reported to contribute to the formation of the iBRB. These proteins include desmocollin-2 (DSC2), protocadherin 12 (PCDH12), fibronectin leucine rich transmembrane protein 3 (FLRT3), latrophillin-3 (ADGRL3/LPHN3), and cell adhesion molecule-1 (CADM1). To validate the model, the presence of previously reported and novel junctional proteins within the retinal microvasculature was confirmed using immunohistochemistry (Figure 3B). It is possible that some of the proteins detected in RMECs are involved in junctions with other cell types, such as pericytes or glial cells in vivo.

The availability of expression data from each cell enables the assessment of not only the average expression of a gene, akin to bulk RNA-seq, but also of the distribution of expression values across the population. Analysis of the variation in expression of genes involved in formation of the microvascular barrier between cells showed that some genes are high in all cells and some low in all cells. Others exhibit a pattern suggestive of bimodal expression, low or undetected in one population of cells and highly expressed in another. Examples of potentially bimodally expressed genes include junctional adhesion molecule 3 (JAM3; Gene ID 83700, OMIM 606871) and Desmocollin 2 (DSC2, Gene ID 1824, OMIM 125645; Figure 4). The modality of each gene can be visualized using violin plots of the expression of the gene across all cells with each dot representing a single cell, as shown in Figure 4A. Dots are sized according to the total reads per cell, and their even distribution demonstrates that normalization has minimized potential technical artifacts due to the variable read number.

In addition to determining expression modality, which provides insights into how genes are regulated, single-cell resolution enables correlation of the expression patterns of different genes between cells. A subset of barrier genes has highly correlated gene expression across individual RMECs (Figure 4B). These genes can be divided into two subgroups within which expression is highly positively correlated but between which there is negative correlation. The observation that a group of genes are expressed together suggests that they may be involved in the same function, although not necessarily physically interacting. A small group of genes centered on PECAM-1 (Gene ID 5175, OMIM 173445) and including CDH5 (Gene ID 1003, OMIM 601120), MCAM (Gene ID 4162, OMIM 155735), ESAM (Gene ID 90952, OMIM 614281), and ADGRL4 (Gene ID 64123, OMIM 616419) stand out for having highly correlated expression. As expected for any subgroup of barrier-related genes, these two sets of coexpressed genes are connected in a network of known protein–protein interactions. Although correlated gene expression is extremely weak evidence for functional interactions, the possibility that these two groups may reflect functional relationships is supported by the clustering of the two sets of proteins in distinct regions of the protein–protein interaction network (Appendix 7).

Cell-specific and novel alternative splicing

In addition to analysis of the gene expression levels between individual cells, scRNA-seq can be used to assess variations in alternative splicing. To achieve this, alternative splicing events within single cells were detected using the Python program Outrigger [31]. The detection of >4,000 skipped exon events and >400 mutually exclusive events emphasizes the diversity of splice isoforms present within RMECs. Although for most events almost all cells have the exon either spliced in (PSI = 1) or spliced out (PSI = 0), there are some for which both isoforms are common, as illustrated by a plot of the global average PSI values for all genes (Figure 5). Differential splicing does not separate the RMECs into obvious clusters, as illustrated with a t-distributed stochastic neighbor embedding (t-SNE) plot and heatmap based on PSI values (Figure 5).

Multiple instances of previously annotated and novel alternative splicing events were observed in the microvascular barrier genes (Appendix 8). For example, the penultimate exon of adducin 1 (ADD1; Gene ID 118, OMIM 102680), which encodes the alpha subunit of a family of cytoskeletal proteins, is skipped in some transcripts that encode a truncated protein with 11 different amino acids at the C-terminus (Figure 6A). This almost certainly has functional consequences due to the loss of phosphorylation sites and the domain involved in interaction with calmodulin [54]. Remarkably, although both isoforms were present in some cells, many expressed one variant exclusively. Similarly, an exon in the adducin 3 (ADD3; Gene ID 120, OMIM 601568) gene was skipped in all transcripts in some cells (i.e., percent spliced in or PSI = 0) and present in all transcripts in other cells (PSI = 1). The functional consequences of the insertion of 32 amino acids toward the C-terminus of the encoded protein due to inclusion of the exon are unknown. Inclusion of an additional exon in transcripts of ECT2 (Gene ID 1894, OMIM 600586), a guanine nucleotide exchange factor (GEF) that participates in the formation of epithelial TJs, also maintains the reading frame and results in the addition of 31 amino acids of unknown function (Figure 6C). In addition to events previously annotated in cattle (or mice or humans), many novel events were observed (for many of which further investigation revealed supporting EST evidence). For example, the inclusion of an exon in the SNAP23 (Gene ID 8773, OMIM 602534) transcript likely has a regulatory function mediated by its extension of the 5′-untranslated region (UTR). Again, the exclusive expression of one SNAP23 isoform in many cells is supported by multiple reads. Reverse transcription (RT)–PCR with primers spanning the exons of interest validated the existence of the alternative splicing events described above (Figure 6E). The PSI values across all cells for the alternative splicing events detected in the BRB model genes can be visualized as a heatmap (Figure 6F). This shows that although most events are either present or absent in most cells, for a significant minority there are both forms in the same cell.