Transcriptomic Responses of Human Retinal Vascular Endothelial Cells to Inflammatory Cytokines.

Purpose: Molecular profiling of human retinal endothelial cells provides opportunities to understand the roles of this cell population in maintenance of the blood-ocular barrier, and its involvements in diverse retinal vasculopathies. We aimed to generate a transcriptome of human retinal endothelial cells in the unstimulated state, and following treatment with inflammatory cytokines linked to cell dysfunction.
Methods: Endothelial cells were isolated from retinae of five human cadaveric donors, and treated for 60 minutes and 24 hours with interleukin-1β or tumor necrosis factor-α, or exposed to medium alone for the same intervals. Expression of intercellular adhesion molecule-1 was measured by RT-qPCR to confirm cytokine-induced activation of the cells. RNA was sequenced on the Illumina NovaSeq 6000 platform. Reads were aligned to the human GRCh38 genome, and reads that aligned to Ensembl-annotated genes were counted. Quality control of sequencing was performed with FastQC, and sequences were classified by Kraken.
Results: A human retinal endothelial cell RNA-sequencing dataset with mean of 99% reads aligned to the human genome was produced as raw RNA sequence data (FASTQ files) and processed read data (XLSX files). Multidimensional scaling analysis showed a strong donor effect, which was readily controlled by ComBat. Conclusions: Our dataset may be useful for human retinal endothelial cell transcriptomic assemblies, functional gene annotating and/or gene expression and enrichment analyses, as well as cross-dataset harmonization. Translational Relevance: The molecular profile of the human retinal endothelium is a source of candidate biologic targets for retinal vasculopathies.

Introduction

Molecular profiling of human cell populations is providing important opportunities to understand normal physiology, dissect the pathophysiology of common and rare diseases, and identify candidate biologic drug targets for those conditions. The endothelial cells that line the vasculature of the retina play an essential role in the blood-ocular barrier, protecting the neural tissue, and facilitating nutrient and metabolite exchange. They also participate in the development of diverse retinal vasculopathies, including diabetic retinopathy, retinal arterial and venous occlusions, retinal infections and non-infectious posterior uveitis. Thus, datasets that describe human retinal endothelial cells under different experimental conditions are of considerable scientific and medical interest.

Multiple technologies exist for investigating the molecular phenotype of a cell, which may include descriptions of the epigenome, various transcriptomes, the proteome and the metabolome. Overall, transcriptomic profiling conducted by RNA-sequencing is the most widely applied approach today, related to considerations that include its discovery focus, high-throughput analysis pipeline, range of sequencing and analysis tools, and relatively low cost, as well as the comprehensive coverage of the cell's molecular status that this advanced technology offers. Relevance for retinal biology and disease have been highlighted by Farkas and colleagues. Already, RNA-sequencing has been used to explore the phenotype of human retinal endothelial cells under different conditions.–

We used RNA-sequencing to profile the transcriptome of endothelial cells isolated from the retinae of five human cadaveric donors; these isolates were processed separately, and individually sequenced. The endothelial cells were studied in the unstimulated state, and following brief and extended treatments with the cytokines, interleukin (IL)-1β or tumor necrosis factor (TNF)-α. These master inflammatory cytokines have been linked to retinal endothelial cell dysfunction in ischemic and inflammatory retinal vasculopathies., The dataset described in this paper is publicly accessible, and its use is not restricted.

Experimental Samples

Five retinal endothelial cell isolates were prepared separately from human cadaver donor eyes, using the positive selection procedure that we have described comprehensively in Methods in Molecular Biology, and used in diverse, published research.– These individual isolates were stored in liquid nitrogen prior to transcriptomic profiling.

The retinal endothelial cells were seeded for confluence in MCDB-131 medium (Merck, Sigma-Aldrich, St. Louis, MO) supplemented with 10% fetal bovine serum (FBS) (Thermo Fisher Scientific-Gibco, Grand Island, NY) and endothelial growth factors (EGM-2 SingleQuots supplement, omitting FBS, hydrocortisone and gentamicin; Clonetics, Lonza, Walkersville, MD) in 12-well plates (Growth area: 4 cm2), rested overnight at 37° C with 5% CO2 in air, and subsequently treated with IL-1β (Bio-Techne-R&D Systems, Minneapolis, MN; Catalogue number: 201-LB-005; Working concentration: 5 ng/mL) or TNF-α (Bio-Techne-R&D Systems; Catalogue number: 210-TA-020; Working concentration: 10 ng/mL) in fresh medium, or fresh medium alone, for 60 minutes or 24 hours. The human cadaver donor eyes used for this work were obtained from the Eye Bank of South Australia (Adelaide, Australia) under an application approved by the Southern Adelaide Clinical Human Research Ethics Committee (Application number: 175.13).

Methods

RNA-Sequencing

Total RNA was recovered from individual endothelial cell cultures using Trizol Reagent (Thermo Fisher Scientific-Invitrogen, Carlsbad, CA), with RNase-free glycogen added during the precipitation to maximize yield, and frozen at −80 °C prior to use for RNA sequencing. An RNA integrity number of at least 8 was confirmed on the LabChip GX Touch 24 Nucleic Acid Analyzer (Perkin Elmer, Hopkinton, MA). The cDNA libraries were prepared using the TruSeq Stranded Total RNA Library Preparation Kit, following depletion of ribosomal RNA with the Ribo-Zero Plus rRNA Depletion Kit (both from Illumina, San Diego, CA). Fragment size was checked on the LabChip GX Touch 24 Nucleic Acid Analyzer, and concentration was determined with the Qubit 2.0 (Thermo Fisher Scientific, Waltham, MA). Equimolar pooled libraries were sequenced on the Illumina NovaSeq 6000 using 2 × 150 paired-end reads, aiming for an average of 50 million reads per sample. The PhiX Control v3 library (Illumina) was used as a sequencing control.

Data Processing

Sequences were trimmed by Trimmomatic v0.38, using the sliding window method, with a window size of 4 and a Phred quality score of 30. Cutadapt v1.17 was used to remove the Illumina universal adapter sequence. Reads were aligned to the GRCh38 human genome with HISAT2 v2.1.0, using strand-specific flag and otherwise default parameters. featureCounts v1.5.0-p2 was used to count strand-specific reads aligned to genome features (genes) based on Ensembl v93 annotation.

Public Repository Link

These RNA sequencing data are stored in the Gene Expression Omnibus (GEO) of the National Center for Biotechnology Information (NCBI). The GEO Series link for this dataset is: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE144785. The GEO recommends linking to the Series record as this provides a summary of the experiment, as well as links to other relevant accessions. The processed read data are provided in Supplementary Table S1.

Data Summary

Donor age and sex, time from death to cell isolation, and cell passage at the time of study, as well as the alignment statistics for each human retinal endothelial cell isolate, are presented in Table 1. The authors and the GEO place no restrictions on the distribution or use of the RNA sequencing data. Published works using these data should cite this publication, the GEO series accession number (GSE144785) and the original GEO database publication or the latest GEO database publication as listed in the GEO.

Table 1.

Description of Human Retinal Endothelial Samples, Including Donor Demographics and Alignment Statistics

Isolate (Passage)	Donor: Sex, Age	Death to Isolation Time	Condition: Cytokine, Time	Paired Reads	Paired Reads Aligning to Genome (%)	Paired Reads Aligning Unambiguously to Annotated Genes (%)
1 (2)	F, 42 yr	20 hr	Control, 60 min	91,880,923	91,228,568 (99.29%)	50,592,086 (55.06%)
			Control, 24 hr	58,665,969	58,425,439 (99.59%)	35,353,899 (60.26%)
			IL-1β, 60 min	54,392,166	54,185,476 (99.62%)	32,791,528 (60.29%)
			IL-1β, 24 hr	61,284,640	60,996,602 (99.53%)	38,127,340 (62.21%)
			TNF-α, 60 min	58,674,151	58,451,189 (99.62%)	40,374,450 (68.81%)
			TNF-α, 24 hr	43,937,185	43,629,625 (99.30%)	24,876,879 (56.62%)
2 (2)	F, 59 yr	10 hr	Control, 60 min	57,203,198	56,831,377 (99.35%)	39,232,816 (68.59%)
			Control, 24 hr	57,061,723	56,662,291 (99.30%)	36,704,318 (64.32%)
			IL-1β, 60 min	67,527,485	67,270,881 (99.62%)	43,837,000 (64.92%)
			IL-1β, 24 hr	54,150,960	53,934,356 (99.60%)	35,183,316 (64.97%)
			TNF-α, 60 min	59,448,395	59,204,657 (99.59%)	34,291,555 (57.68%)
			TNF-α, 24 hr	44,707,449	43,706,002 (97.76%)	28,200,494 (63.08%)
3 (3)	F, 43 yr	24 hr	Control, 60 min	57,916,040	57,695,959 (99.62%)	38,687,922 (66.80%)
			Control, 24 hr	50,261,528	50,045,403 (99.57%)	33,012,471 (65.68%)
			IL-1β, 60 min	66,313,255	65,968,426 (99.48%)	41,168,092 (62.08%)
			IL-1β, 24 hr	50,499,562	50,247,064 (99.50%)	32,402,699 (64.16%)
			TNF-α, 60 min	60,050,927	59,978,866 (99.88%)	9,897,383 (16.48%)
			TNF-α, 24 hr	50,292,931	50,031,408 (99.48%)	31,064,967 (61.77%)
4 (2)	F, 54 yr	31 hr	Control, 60 min	55,981,156	55,757,231 (99.60%)	37,892,323 (67.69%)
			Control, 24 hr	54,062,301	53,802,802 (99.52%)	33,214,431 (61.44%)
			IL-1β, 60 min	63,682,704	63,427,973 (99.60%)	34,960,173 (54.90%)
			IL-1β, 24 hr	48,600,754	48,391,771 (99.57%)	31,419,265 (64.65%)
			TNF-α, 60 min	60,328,580	59,791,656 (99.11%)	33,880,535 (56.16%)
			TNF-α, 24 hr	72,695,997	72,019,924 (99.07%)	46,217,221 (63.58%)
5 (2)	F, 59 yr	18 hr	Control, 60 min	51,457,569	51,236,301 (99.57%)	34,115,169 (66.30%)
			Control, 24 hr	48,618,962	48,385,591 (99.52%)	33,696,362 (69.31%)
			IL-1β, 60 min	60,230,290	60,079,714 (99.75%)	13,097,586 (21.75%)
			IL-1β, 24 hr	49,383,724	49,126,929 (99.48%)	32,349,949 (65.51%)
			TNF-α, 60 min	55,171,735	54,796,567 (99.32%)	35,131,641 (63.68%)
			TNF-α, 24 hr	47,890,187	47,636,369 (99.47%)	30,255,012 (63.18%)

Abbreviations: F = female, IL-1β = interleukin-1β, TNF-α = tumor necrosis factor-α, yr = years, min = minutes, hr = hours.

Technical Validation

Endothelial Activation

Endothelial cell activation is characterized by the induction of adhesion molecules, including the immunoglobulin superfamily member, intercellular adhesion molecule (ICAM)-1. Prior to performing the RNA sequencing, ICAM-1 transcript was measured on total cellular RNA by reverse transcription (RT) and real-time quantitative polymerase chain reaction (qPCR), to confirm activation of the endothelium by IL-1β and TNF-α for both treatment intervals.

Total RNA was reverse transcribed using iScript Reverse Transcription Supermix for RT-qPCR (Bio-Rad, Hercules, CA). Total RNA input was 125 ng per reaction, yielding 20 µl of cDNA which was diluted 10-fold ahead of use in qPCR. Quantitative real-time PCR was performed using the CFX Connect Real Time PCR Detection System (Bio-Rad), and a reaction mix consisting of 4 µL of iQ Sybr Green Supermix or SsoAdvanced SYBR Green Supermix (Bio-Rad), 1.5 µL each of 20 µM forward and reverse primer, 11 µL of nuclease-free water and 2 µL of diluted cDNA. Amplification consisted of: a pre-cycling hold at 95 °C for 5 minutes; 40 cycles of denaturation for 30 seconds at 95 °C; annealing for 30 seconds at 60 °C; and extension for 30 seconds at 72 °C. Relative expression of ICAM-1 was determined by the Pfaffl method, and normalized to ribosomal protein lateral stalk subunit P0 (RPLP0), which was stable across treatments and time points according to the coefficient of variance.

Primer sequences (and expected product sizes and efficiencies) were as follows:

ICAM-1, forward 5’-TAAGCCAAGAGGAAGGAGCA-3’,

reverse 5’-CATATCATCAAGGGTTGGGG-3’ (289 bp, 100.8%);

and RPLP0, forward 5’-GCAGCATCTACAACCCTGAA-3’,

reverse 5’-GCAGATGGATCAGCCAAGAA-3’ (235 bp, 91.8%).

This assessment of the human retinal endothelial cell isolates by RT-qPCR demonstrated a statistically significant increase in expression of ICAM-1 transcript following treatment with both IL-1β and TNF-α for 60 minutes and 24 hours, with greater differences at 24 hours, confirming effective activation of the endothelial cells (Figs. 1A and 1B). The cytokine-induced increase in ICAM-1 expression in the cells was measured at a similar level, and with identical significance, by RNA-sequencing (Figs. 1C and 1D).

Figure 1.

Technical validations of the transcriptomic dataset generated for human retinal endothelial cells treated with interleukin (IL)-1β or tumor necrosis factor (TNF)-α. (A and B). Graphs showing relative transcript expression for intercellular adhesion molecule (ICAM)-1 in cytokine-treated cells after (A) 60 minutes and (B) 24 hours in comparison to medium control by RT-qPCR. Reference gene is ribosomal protein lateral stalk subunit P0 (RPLP0). (C and D). Graphs showing ICAM-1 transcript expression by RNA-sequencing for the same samples. Each circle represents one individual isolate, with crossbars indicating mean. n = 5 isolates /condition. Data were analyzed by Friedman test. (E). Bar plot representing RNA composition of medium control and cytokine-treated cells, as classified by Kraken against the NCBI RefSeq database. (F and G). Multidimensional scaling plots showing global gene expression in medium control and cytokine-treated cells (F) pre- and (G) post-adjustment for donor effect with ComBat in svaseq. Colors indicate individual isolates; shapes indicate treatments; open shapes indicate 60-minute treatments and closed shapes indicate 24-hour treatments; circles represent groupings based on individual isolates. Additional abbreviations: C = control; T = TNF-α; I = IL-1β; PC = principal component; min = minutes; hr = hours.

RNA Quality

Prior to further processing, quality control of the sequenced RNA was performed with FastQC v0.11.3 and summarized with MultiQC v1.7. The number of reads per sample averaged 56.8 million. Mean Phred scores across all bases for both forward and reverse reads exceeded 35 (i.e. probability of incorrect bases approximated 1 in 3200), which confirmed that the sequence data were of high quality. Mean ribosomal RNA content in the 30 samples was estimated to be 10.1% by SortMeRNA v2.1b, indicating effective depletion.

Mycoplasma Contamination

Contamination of cultured mammalian cells with Mycoplasma species occurs commonly, and the presence of these bacteria broadly impacts cell biology, including gene expression changes, with potential to produce misleading results. Prior to alignment, RNA sequences were classified against the NCBI RefSeq database using Kraken v2.0.7. Assessment of the 30 samples demonstrated that a mean of 99.4% of reads aligned to the human genome, and no reads aligned to Mycoplasma sequence at the 1% classification limit of the software tool (Fig. 1E).

Biological Variability

The EdgeR v3.26.5 package from Bioconductor was used to normalize read counts for library size with the trimmed mean of M-values (TMM) method, and to perform multidimensional scaling analysis. This analysis revealed a donor effect on gene expression (Fig. 1F). This is expected in any study involving samples taken from multiple human donors,– and was observed in a previous oligonucleotide microarray study of the human retinal endothelial cell transcriptome. However, as our experimental design included both treated and untreated samples for each isolate, this donor effect could be readily controlled by ComBat implemented in the svaseq Bioconductor package v3.38.0 (Fig. 1G).

Discussion

Our report describes a publicly available dataset of the transcriptome of human retinal endothelial cells in resting state, and IL-1β- and TNF-α-activated states, generated across five human donors. These data provide the opportunity for differential expression analysis – with follow-up Gene Ontology, pathway and network enrichment analyses – to study the molecular phenotype of the human retinal endothelial cell and its response to two master inflammatory cytokines. Additional enrichment analyses could explore regulatory processes, focusing on discovery of highly active signaling pathways, identification of connected long noncoding RNAs and prediction of central transcription factors. Multiple tools are publicly available for these purposes, e.g. InnateDB, ChEA3, LncSEA, and SPAGI.

Given the roles of IL-1β and TNF-α in retinal endothelial dysfunction, analyses using our dataset may be used to illuminate basic mechanisms of retinal vasculopathies, ranging from diabetic retinopathy to posterior uveitis. Apart from providing the comparison with cytokine-treated cells, control samples provide data that separately build understanding of the unique molecular phenotype of the human retinal endothelial cell, including its involvement in the blood-ocular barrier. Research using our dataset may be expanded by comparisons with relevant datasets that have been provided publicly by ourselves and others.

As molecular profiling technologies have advanced, we have generated a series of datasets using human retinal endothelial cells, which have been shared as supplements to publications and/or deposits in public repositories that are identified in the publications. These datasets include: (1) Affymetrix Human Genome Focus Array oligonucleotide microarray profiles of human retinal endothelial cells in comparison to donor-matched choroidal endothelial cells, stimulated with lipopolysaccharide or Toxoplasma gondii; (2) comparison of protein expression differences between human retinal and choroidal endothelial cells by two-dimensional difference gel electrophoresis tandem mass spectrometry, and (3) shotgun reverse-phase liquid chromatography tandem mass spectrometry profile of human retinal endothelial cells in isolation and in comparison to donor-matched choroidal endothelial cells., We have also employed Bio-Rad PrimePCR Reference Genes H96 and H96Plus PCR arrays to identify stable reference genes for use in RT-qPCR studies of human retinal endothelial cells.

Since our human retinal endothelial cell isolation, culture and stimulation techniques are standard, the datasets that we have produced provide opportunity for cross-dataset harmonization, for both biologically focused research and to compare outputs from different profiling methodologies. However, there are caveats to such comparisons., Microarray detection is limited by the probe series, which in our study detected approximately 8500 transcripts, compared to over 11,000 transcripts discovered by RNA-sequencing in this work. Additionally, overall microarray has a smaller dynamic range and thus may be less accurate for quantifying gene expression levels. Abundance of transcript and protein may not correlate well, due to biological factors including turn-over, and translation-related regulation, as well as technical factors. However, we have reported good agreement between the human retinal endothelial transcriptome and proteome by focusing on molecular groupings, rather than individual molecules.

Several other groups have used RNA-sequencing to generate different datasets with human retinal endothelial cells. Savage et al., produced two datasets based on these cells, comparing the effect of transcription factor blocking agents, INCA-6 or GSK0660, during a 4-hour exposure to TNF-α, as well as no cytokine treatment. Robinson et al. reported a dataset for cells activated by overnight treatment with IL-6 and soluble IL-6 receptor, versus controls that were unmanipulated or exposed to the IL-6 signaling blocker, sgp130Fc. Huang et al. described a dataset from cells treated for 48 hours with placental growth factor antibody versus non-treated control cells. Datasets for these studies were deposited in the Sequence Read Archive hosted by NCBI. Shao et al. generated an RNA-sequencing dataset for human retinal endothelial cells treated with different glucose concentrations and transthyretin for 48 hours, and deposited this in GEO. One caveat in comparing these datasets with the one we report is that the other groups sourced the human retinal endothelial cells commercially – or in the study by Shao et al., the source was not identified – and used them en masse. Our dataset includes separately sequenced and processed data for five different human donors. Browning et al. also performed a multi-donor study, comparing human retinal, choroidal and iris endothelial cells from three donors, as well as human umbilical vein cells; however, in this work, profiling was by oligonucleotide microarray.

The dataset reported by Shao et al. was presented in a format that allowed ready extraction of read counts. Representing transcript levels relative to RPLP0, which we have previously identified as a stable reference gene for expression studies using human retinal endothelial cells, permitted a direct comparison with our dataset. This comparison is presented in Table 2. Gene expression of untreated, control cells showed similar expression across the two datasets for 13 of 15 common endothelial cell transcripts, including vascular endothelial (VE)-cadherin, cell adhesion molecules from the immunoglobulin superfamily, and receptors for growth factors such as vascular endothelial growth factor. Two transcripts had obviously higher expression in Shao et al's dataset (activated leukocyte cell adhesion molecule [ALCAM] and epidermal growth factor [EGFR]), for which there are numerous potential explanations, including the use of a single donor in that study (versus five in ours), and bench and bioinformatic methodological differences, including those related to cell isolation and culture, RNA extraction, library preparation and sequencing, read alignment and counting, and/or versions of algorithms and databases.

Table 2.

Expression of Selected Retinal Endothelial Cell Transcripts Under Control (no Treatment) Conditions in This Study and a Comparable RNA-Sequencing Study. Transcripts Are Identified by Gene Name, and Expression is Presented Relative to Ribosomal Protein Lateral Stalk Subunit P0 (RPLP0, ENSG00000089157)

	Mean Normalized Expression
Transcript*	Ryan et al., 2022 (n = 5 Isolates)	Shao et al., 20198 (n = 3 Replicates)
VWF	1.143915	1.800747
PECAM1	1.118422	1.620269
ICAM1	0.747052	0.275639
VCAM1	0.5047	0.483821
ALCAM	0.850183	6.058259
SELE	0.504707	0.660358
CDH5	1.035694	0.736351
FLT1	0.777363	1.506564
KDR	0.968024	1.161705
FLT4	0.604592	1.149256
TIE1	0.857932	0.772019
TEK	0.713154	0.754045
IGF1R	0.52761	0.622734
FGFR1	0.620695	1.563234
EGFR	0.484111	12.52205

Abbreviations: VWF = von Willebrand factor, PECAM1 = platelet and endothelial cell adhesion molecule 1, ICAM1 = intercellular adhesion molecule 1, VCAM1 = vascular cell adhesion molecule 1, ALCAM = activated leukocyte cell adhesion molecule, SELE = selectin E, CDH5 = cadherin 5, FLT1 = fms related receptor tyrosine kinase°1, KDR = kinase insert domain receptor, FLT4 = fms related receptor tyrosine kinase 4, TIE1 = tyrosine kinase with immunoglobulin like and EGF like domains°1, TEK = TEK receptor tyrosine kinase, IGF1R = insulin like growth factor 1 receptor, FGFR1 = fibroblast growth factor receptor 1, EGFR = epidermal growth factor receptor.

While studies using human cells require some degree of in vitro manipulation, particularly when different conditions will be compared, one limitation across all such work is the potential for cell phenotype to be impacted. An addition caveat to consider in all the studies with human retinal endothelial cells, including ours, is that the primary cell isolates include a mixed population of human arterial, venous and capillary retinal endothelial cells. In using this dataset to progress research on the human retinal endothelium, these limitations can be balanced by validating leads across the independent, open access datasets, as well as in follow-up studies with additional cell isolates and/or appropriate in vivo experimental animal models.

Data Repository

These data are stored in text-based open formats and do not require a dictionary to read: (1) raw RNA sequence data in FASTQ file format; and (2) processed read data in XLSX file format. The files may be directly downloaded for use in transcriptome assemblies or gene expression analyses. There are no missing values in these files.