Molecular Vision 2020; 26:766-779
Received 17 February 2020 | Accepted 25 November 2020 | Published 27 November 2020
1Wellcome-Wolfson Institute for Experimental Medicine, Queen’s University Belfast, Belfast, Northern Ireland; 2Genomics Core Technology Unit, Faculty of Medicine, Health & Life Sciences, Queen’s University Belfast, Belfast, Northern Ireland
Correspondence to: David Simpson, The Wellcome – Wolfson Institute for Experimental Medicine, School of Medicine, Dentistry and Biomedical Sciences, Queen’s University Belfast, 97 Lisburn Road, Belfast BT9 7BL. Phone: +44-28- 90976470; email: David.Simpson@qub.ac.uk
Purpose: To better characterize retinal endothelial barrier properties through analysis of individual transcriptomes of primary bovine retinal microvascular endothelial cells (RMECs).
Methods: Individual RMECs were captured on the Fluidigm C1 system, cDNA libraries were prepared using a Nextera XT kit, and sequencing was performed on a NextSeq system (Illumina). Data analysis was performed using R packages Scater, SC3, and Seurat, and the browser application Automated Single-cell Analysis Pipeline (ASAP). Alternative splicing events in single cells were quantified with Outrigger. Cytoscape was used for network analyses.
Results: Application of a single-cell RNA sequencing (scRNA-seq) analysis workflow showed that RMECs form a relatively homogeneous population in culture, with the main differences related to proliferation status. Expression of markers from along the arteriovenous tree suggested that most cells originated from capillaries. Average gene expression levels across all cells were used to develop an in silico model of the inner blood–retina barrier incorporating junctional proteins not previously reported within the retinal vasculature. Correlation of barrier gene expression among individual cells revealed a subgroup of genes highly correlated with PECAM-1 at the center of the correlation network. Numerous alternative splicing events involving exons within microvascular barrier genes were observed, and in many cases, individual cells expressed one isoform exclusively.
Conclusions: We optimized a workflow for single-cell transcriptomics in primary RMECs. The results provide fundamental insights into the genes involved in formation of the retinal–microvascular barrier.
The inner blood–retina barrier (iBRB) controls fluid exchange across retinal capillary beds. Tight junctions (TJs) and adherens junctions (AJs) between adjacent vascular endothelial cells result in high transendothelial electrical resistance (TEER) and greatly restrict the paracellular flow of water, proteins, lipids, and immune cells into the retina . Breakdown of the iBRB can have serious pathophysiological consequences in the retina and causes the edema that contributes to vision loss and blindness in diabetic retinopathy, retinal vein occlusions, retinopathy of prematurity, and uveitis .
The TJ proteins that make up the iBRB include members of the tetraspanin (claudins, occludin, tricellulin, and marvelD3 proteins) and junctional adhesion molecule (JAM) families . These proteins interact with cytosolic scaffolding proteins, such as zonula occludens (ZO-1), afadin (AF6), and cingulin, which, in turn, anchor the TJ complex to the actin cytoskeleton . The main AJ protein of the iBRB is vascular endothelial (VE)-cadherin, which is a member of the classical cadherin superfamily . VE-cadherin complexes with β-catenin, which links through other proteins to the actin cytosketeton . However, the specific barrier genes employed by retinal endothelial cells, their interactions with one another, and the heterogeneity in their expression between cells are not well understood. The recent development of single-cell transcriptomic approaches now provides an opportunity to investigate these factors, all of which have consequences for the development of strategies for regulating permeability.
Traditional bulk RNA sequencing or microarray approaches compare average gene expression between populations of cells and are unable to discriminate heterogeneity among individual cells. In contrast, single-cell RNA sequencing (scRNA-seq) technology [6-10] measures the gene expression profiles of individual cells and has enabled characterization of existing and new retinal cell types at the transcriptome level in adult  and developing retinas . The cellular resolution provided by scRNA-seq enables coexpression analysis between cells to suggest potential functional interactions between genes.
In many cases, the function of a gene depends not only on its level of expression but also on the specific alternatively spliced variants present, which can give rise to different protein isoforms. Alternative splicing has been recognized as a regulatory force in TJ assembly for almost two decades . Several ZO-1 splice isoforms have been localized to TJs , and the ZO-1 alpha(+) variant has been shown to enhance tightness, while the ZO-1 alpha(−) variant is present in dynamic junctions that are readily opened by physiologic signals [14,15]. Occludin 1B’s wide epithelial distribution and conservation across species suggest a potentially important role in the structure and function of the TJ . An alternatively spliced occludin isoform generated by skipping exon 4 has been identified , and four differentially spliced occludin mRNA transcripts have been reported in human epithelial tissues (the placenta and the colon) and colonic epithelial cells . MarvelD3 colocalizes with occludin at TJs  and is expressed as two alternatively spliced isoforms by different types of epithelial as well as endothelial cells. TJ protein claudin-10 exists in two isoforms with alternative exons, 1a and 1b (Cldn10a, Cldn10b) with Cldn10a restricted to the kidney . Four additional claudin-10 splice variants have since been found in mice and humans . Some isoforms insert into the TJ, whereas those lacking exon 4 are retained in the endoplasmic reticulum.
Bulk RNA sequencing has demonstrated the rich diversity of alternatively spliced transcript variants [22,23], but it cannot determine whether cell populations are homogeneous or heterogeneous . This knowledge is important to identify previously unrecognized subpopulations expressing distinct splice isoforms and to determine the extent of differential isoform expression that might change in response to environmental conditions. For example, are there subpopulations of endothelial cells expressing different isoforms of junction proteins, and does differential isoform expression support alternative splicing as a possible mechanism for regulating cell junction properties? Single-cell RNA sequencing has revealed heterogeneity of alternative splicing between the individual cells of many populations [25-30], but it is unclear to what extent all cells in a population express only one transcript variant, one of several variants (bimodality), or simultaneous expression of multiple transcript variants.
Although TJ and AJ transcripts can undergo alternative splicing giving rise to multiple protein isoforms, there have been no studies specifically investigating the splicing patterns of iBRB genes. Single-cell RNA-seq can be used to assess heterogeneity in overall gene expression and the choice of splice isoforms between individual cells .
We report the application of scRNA-seq to characterize the transcriptomes of primary bovine retinal microvascular endothelial cells (RMECs). It would be difficult to isolate a pure population of RMECS directly from retinal tissue, and primary, early passage cultured RMECs provide a widely used model for examining the iBRB that provides the practical utility of working with cultured cells while maintaining as many in vivo features of the cells as possible. These data underpin a retina-specific model of endothelial barrier function, including widespread cell-specific patterns of alternative splicing.
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 . 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 cm2, 1% gelatin-coated culture flask (TPP, Trasadingen, Switzerland) and maintained at 37 °C in 5% CO2. 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 Na2HPO4, 1.8 mM KH2PO4, 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.
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 CO2 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 Na2HPO4, 1.5 mM KH2PO4) 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 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.
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 × 105 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 was performed by the Genomics Core Technology Unit (GCTU) unit at Queen’s University Belfast (Genomics) using the NextSeq® 500 system (Illumina).
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 , and cells with <100,000 aligned reads were discarded. Alternative splicing within single cells was quantified using Outrigger , 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)  was used to visualize splicing and create sashimi plots of alternative splicing events. Gene counts were calculated using HTseq , and an expression matrix was created. Quality control (QC) was visualized at each stage using multiQC . Seurat  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)  was used to generate violin plots of gene expression. The single-cell latent variable model (scLVM)  within the web-based single-cell transcriptomic analytical tool Automated Single Cell Analysis Pipeline (ASAP)  was used for normalization of gene expression, accounting for differences in sequencing depth and technical bias. Single-cell consensus clustering (SC3)  was used for unsupervised clustering of the filtered data set within ASAP . 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 , 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  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  was used in conjunction with the GENEmania  and String  applications for network analyses and generation of the blood–retina barrier model. Protein functions and interactions were analyzed using the STRING database . CorrelationCalculator was used to calculate and generate heatmaps of pairwise correlations between genes based on their expression pattern among individual cells . Correlation results were imported into Cytoscape and MetScape used to view correlation networks [51,52].
In this study, we implemented an optimized workflow for scRNA-seq analysis of RMECs (Appendix 3). Two replicate batches of cells were processed, with 49 and 51 single cells, respectively, passing QC and taken through the downstream analysis pipeline. The raw FastQ files and expression matrices are available from Gene Expression Omnibus (GEO) accession number GSE136800.
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  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  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 . 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).
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).
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 . 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 . 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.
To the best of our knowledge, this is the first investigation of primary retinal endothelial gene expression at the level of individual cells. Although optimized for RMECs, the experimental and bioinformatics approach we developed would be equally applicable to cultures of other retinal cell types. Comparison of individual RMEC transcriptomes demonstrated that they form a relatively homogeneous population in culture, with the main variation related to the proliferative state rather than the genotype. The regressing out of these cell cycle effects enabled more cells to be retained within subsequent analyses. Expression of in vivo markers suggests that the RMECs most closely resemble capillary endothelial cells, which may reflect either their origin or convergence on a common phenotype in culture .
We exploited the single-cell transcriptome data to better understand the molecular properties of the retinal endothelial cell barrier. We have made the raw data publicly available (GEO ID GSE136800) so that they can be used by other researchers to focus on other features of retinal endothelial cell biology, such as cell matrix adhesion, receptor signaling pathways, proliferation, and angiogenesis. The combination of single-cell data with known protein–protein interactions enabled us to develop the most comprehensive in silico model of the iBRB published to date. Although the transcriptome data were obtained from a monolayer of RMECs in culture, the predicted expression of key proteins was validated in retinal vessels in vivo. In addition to confirming the involvement of known genes, we identified novel genes not previously associated with the iBRB. The model incorporates previously characterized interactions between all of these junctional and structural proteins including links to the cytoskeleton. The scRNA-seq data enabled us to investigate how these barrier genes were coexpressed between individual cells and identify networks of genes with correlated expression. These results suggest that individual cells differentially regulate specific subgroups of genes that influence their barrier properties. Such an insight would not have been possible without the use of a single-cell analysis approach.
Alternative splicing is a critical mechanism for generating diverse protein isoforms with differing functions. For example, the importance of VEGFA splice isoforms in regulating retinal angiogenesis is well-known , but analysis of alternative splicing of genes involved in the formation of the iBRB has been limited. To date, alternative splicing has been determined mainly from bulk RNA-seq data, but scRNA-seq provides several important advantages. First, it makes it possible to determine whether splice events occur with similar frequency in all cells or more or less frequently in subpopulations of cells. This makes it possible to address the question of whether two splice variants are present together in all cells or are mutually exclusive. Second, rare splice events that would be dismissed as background in bulk data may be revealed by scRNA-seq to be common in a small number of individual cells. Therefore, we exploited these RMEC scRNA-seq data to provide a comprehensive analysis of splicing within the genes in the iBRB model. In addition to known splice events, we identified numerous novel splicing events in barrier genes that have not been previously reported in the literature. Perhaps surprisingly, many events were detected as the dominant isoform in a subset of cells while being absent in others, although this is in agreement with previous reports of the dominance of a single transcript variant in each cell . However, we did not detect subpopulations of RMECs based on differential expression of alternative splice isoforms. Nonetheless, these data highlight the heterogeneity in splicing of barrier genes at the cell level that must be considered if we wish to fully understand the molecular properties of the iBRB.
In conclusion, we developed and optimized workflows for studying single-cell transcriptomics in primary RMECs. The gene expression and splicing data from these cells have been made publicly available for secondary analysis. This analysis of barrier genes enabled development of an in silico model of the iBRB, including heterogeneity of gene expression and splicing. This work has enabled a better understanding of the molecular basis of the iBRB and provides an important platform to determine more precisely why this structure breaks down in multiple retinal disorders.
Appendix 1. Primary antibodies used for immunofluorescence.
Appendix 2. Primers used for RT–PCR.
Appendix 3. Single Cell RNA-Seq workflow and analysis pipeline.
Appendix 4. Unsupervised clustering of individual RMECs according to gene expression without cell-cycle correction. (A)
Appendix 5. Genes differentially expressed between proliferating and non-proliferating cells.
Appendix 6. List of barrier genes with alternative IDs and expression values.
Appendix 8. Alternative splicing events within barrier genes.
Single cell libraries were prepared and sequenced with assistance from the Queens University Belfast Genomics CTU. Thanks to Michael O'Hare for assistance with imaging. This work was supported by a Northern Ireland Department for the Economy (DfE) studentship, the Health & Social Care R&D Division, Northern Ireland (STL/4748/13), the Medical Research Council (MC_PC_15026), the Biotechnology and Biological Sciences Research Council (BB/I026359/1), Fight for Sight project grant 5091/5092 and US-Ireland Partnership 14/US/B3116. Sequence data has been deposited in the Gene Expression Omnibus (GEO) repository, accession number GSE136800.