October 2019
Volume 60, Issue 13
Open Access
Retina  |   October 2019
Induced Expression of VEGFC, ANGPT, and EFNB2 and Their Receptors Characterizes Neovascularization in Proliferative Diabetic Retinopathy
Author Affiliations & Notes
  • Yaping Li
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
    Department of Ophthalmology, Second Hospital, Jilin University, Changchun, China
  • Dong Chen
    Center for Genome Analysis, ABLife, Inc., Wuhan, People's Republic of China
    Laboratory for Genome Regulation and Human Health, ABLife, Inc., Wuhan, People's Republic of China
  • Luguo Sun
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Yannan Wu
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Ying Zou
    Department of Ophthalmology, Second Hospital, Jilin University, Changchun, China
  • Chen Liang
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Yongli Bao
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Jingwen Yi
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Yu Zhang
    Center for Genome Analysis, ABLife, Inc., Wuhan, People's Republic of China
  • Jing Hou
    Center for Genome Analysis, ABLife, Inc., Wuhan, People's Republic of China
  • Zhen Li
    Center for Genome Analysis, ABLife, Inc., Wuhan, People's Republic of China
  • Fengyun Yu
    Laboratory for Genome Regulation and Human Health, ABLife, Inc., Wuhan, People's Republic of China
  • Yanxin Huang
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Chunlei Yu
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Lei Liu
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Zaoxia Liu
    Department of Ophthalmology, Second Hospital, Jilin University, Changchun, China
  • Yi Zhang
    Center for Genome Analysis, ABLife, Inc., Wuhan, People's Republic of China
    Laboratory for Genome Regulation and Human Health, ABLife, Inc., Wuhan, People's Republic of China
  • Yuxin Li
    National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, Changchun, China
  • Correspondence: Yuxin Li, National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, 2555 Jingyue Avenue, Changchun 130117, China; liyx486@nenu.edu.cn
  • Yi Zhang, Center for Genome Analysis, ABLife, Inc., 388 Gaoxin 2nd Road, East Lake High-Tech Development Zone, Wuhan 430075, People's Republic of China; yizhang@ablife.cc
  • Zaoxia Liu, Department of Ophthalmology, Second Hospital, Jilin University, Changchun 130021, China; Liuliu15898@sina.com
  • Lei Liu, National Engineering Laboratory for Druggable Gene and Protein Screening, Northeast Normal University, 2555 Jingyue Avenue, Changchun 130117, China; liul905@nenu.edu.cn
Investigative Ophthalmology & Visual Science October 2019, Vol.60, 4084-4096. doi:https://doi.org/10.1167/iovs.19-26767
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Yaping Li, Dong Chen, Luguo Sun, Yannan Wu, Ying Zou, Chen Liang, Yongli Bao, Jingwen Yi, Yu Zhang, Jing Hou, Zhen Li, Fengyun Yu, Yanxin Huang, Chunlei Yu, Lei Liu, Zaoxia Liu, Yi Zhang, Yuxin Li; Induced Expression of VEGFC, ANGPT, and EFNB2 and Their Receptors Characterizes Neovascularization in Proliferative Diabetic Retinopathy. Invest. Ophthalmol. Vis. Sci. 2019;60(13):4084-4096. doi: https://doi.org/10.1167/iovs.19-26767.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
Abstract

Purpose: To investigate whole transcriptional differences between proliferative diabetic retinopathy (PDR) neovascular membranes (NVMs) and retinas, and the regulatory genes participating in retinal neovascularization in PDR.

Methods: We used high-throughput sequencing technology to capture the whole-genome gene expression levels of all participants, including 23 patients with PDR or branch retinal vein occlusion (BRVO), 3 normal retinal samples, and 2 retinal samples from type II diabetic (T2D) eyes by donation, followed by analyses of expression patterns using bioinformatics methods, then validation of the data by in situ hybridization and Western blotting.

Results: We showed that transcriptional profiles of the NVMs were distinct from those of the retinas. Angiogenesis growth factors VEGFC, ANGPT1, ANGPT2, and EFNB2, and their receptors FLT4, TIE1, TIE2, and EPHB4, respectively, were overexpressed. Expression of VEGFA was highly upregulated in T2D retina, but low in the NVMs, while angiogenesis transcription factors, including ETS1 and ERG, were coordinately upregulated in NVMs.

Conclusions: This study described a PDR neovascularization model in which pathological retina-secreted vascular endothelial growth factor A (VEGFA) enhanced the expression of a set of angiogenesis transcription factors and growth factors, to cooperatively induce the retinal neovascularization. Based on these results, novel potential therapeutic targets and biomarkers for PDR treatment and diagnosis are suggested.

Diabetes mellitus is mainly caused by impaired insulin production or action, and does significant long-term damage to many organs, including eyes.1,2 Diabetics have high risks of macrovascular and microvascular complications induced by oxidative stress.3 Diabetic retinopathy (DR) is one of the most common microvascular complications of diabetes.4,5 The advanced stage of DR, proliferative diabetic retinopathy (PDR) with retinal neovascularization, is closely related to vision loss.68 Investigations of the molecular mechanisms involved in retinal neovascularization have shown that vascular endothelial growth factor A (VEGFA), also known as VEGF, is a master causative factor for retinal neovascularization.4,9 VEGF was identified as a diffusible mitogenic heparin-binding protein in the 1980s,10,11 belonging to a gene family that includes VEGFB, VEGFC, VEGFD, and placental growth factor (PlGF).12 These proteins, as well as other vascular endothelium-specific growth factors, play differential and cooperative roles in promoting and regulating angiogenesis or neovascularization, by interacting with their receptors and then stimulating the following processes.1217 
During the pathogenesis of PDR, it has been widely reported that VEGFA-VEGFR2 (VEGF receptor 2, also named KDR or Flk-1) mediates endothelial cell mitogenesis and vascular permeability, via VEGFR2-mediated signaling cascades.9,16 VEGFA could be increased by altered glucose levels, insulin, and other growth factors.1820 The secreted VEGFA level from retinal pigment epithelial cells was increased in the ocular fluid of DR patients and in the retinal ischemic monkey model,2124 and its expression level in DR patients has been correlated with disease progression.24 Consequently, anti-VEGFA therapy has been developed, with the first inhibitor approved in 2004.9,22 The ocular anti-VEGFA therapies have shown clear efficacies in clinical trials.25,26 In addition to VEGFA, other proteins of the VEGF family also have important roles in angiogenesis. Competitive and collaborative relationships among VEGFs and their receptors constitute an important regulatory network during angiogenesis.2731 In human retinal pigment epithelial cells, VEGFA promotes the expression of VEGFC and VEGFR3.31 Therefore, the roles of VEGF family members and angiopoietin and ephrin family members,16,32 as well as their potential cooperation during retinal neovascularization in humans, should be further studied with advanced technologies.33 
An omics study on gene expression alterations of fibrovascular membranes (FVMs) between PDR patients and normal retinal controls has revealed 91 genes that are preferentially expressed in FVMs, when compared with normal retinas.34 In the present study, we obtained qualified cDNA libraries from the small neovascular membranes (NVMs) of patients and retinas from normal or type II diabetic (T2D) individuals, and then sequenced the qualified cDNA libraries (RNA-Seq). The transcriptome profiles of the NVMs were distinct from retinal membranes. We also discovered that a number of genes belonging to the master angiogenesis transcription factor (TF) ETS1 were systematically upregulated during retinal neovascularization. This expression network suggested distinct VEGF regulatory programs, involving VEGFA acting as a trigger of the program switch. 
Materials and Methods
Study Approval, Sample Collection, Treatment, and RNA Extraction
This study was approved by the Ethics Committee of the Second Hospital of Jilin University, and followed the tenets of the Declaration of Helsinki. Informed consent was obtained from the participants after explanation of the nature and possible consequences of the study. NVMs were surgically obtained from 23 eyes of patients with PDR or branch retinal vein occlusion (BRVO) during vitrectomies, which were performed at the Second Hospital of Jilin University. Three normal and two T2D retinas were from donated eyes. Surgical indications included persistent vitreous hemorrhage, NVM formation, and NVM traction retinal detachment. The surgical specimen for each case was frozen in liquid nitrogen immediately after removal. 
After homogenizing the sample with TRIzol reagent, the samples were incubated at room temperature to permit complete dissociation of the nucleoprotein complex. The homogenate was then centrifuged at 4°C, and the supernatant was collected and added to chloroform. The solution was shaken vigorously and centrifuged at 4°C to produce a clear upper aqueous layer (containing RNA), which was precipitated from the aqueous layer with isopropanol. The precipitated RNA was washed with 75% ethanol to remove impurities, and then resuspended with diethyl pyrocarbonate (DEPC)-treated water for use in downstream applications. 
RNA-Seq Library Construction
The quality, integrity, and quantity of the purified RNAs were determined using a 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA). For each sample, 50 ng total RNA was used for the library preparation. Ribosomal RNAs were depleted with a Ribo-Zero rRNA depletion kit for humans (Epicentre, San Diego, CA, USA) before being used for directional RNA-Seq library preparation. Purified RNA was fragmented, and then reverse transcribed with random primers containing adaptor sequences and a randomized hexamer. The cDNAs were then synthesized with terminal-tagging oligo-cDNA using the ScriptSeq v2 RNA-Seq Library Preparation Kit (Epicentre). The cDNAs were purified and amplified. PCR products corresponding to 300 to 500 base pair (bp) were purified, quantified, and stored at −80°C until sequenced. 
The libraries were prepared following the manufacturer's instructions and applied to an Illumina NextSeq 500 platform for 151 nucleotide pair-end sequencing (ABLife, Inc., Wuhan, China). The RNA-Seq data have been deposited in the GEO database with accession number GSE102485. 
Data Processing and Differentially Expressed Gene (DEG) Analysis
Cutadapt35 was used to trim adaptors (with a 10% error) and low-quality bases (less than 16) of raw reads. Reads less than 16 nucleotides in both ends were discarded. We used TopHat236 software to align the filtered reads to the human reference genome (GRCH38) with four read mismatches using an end-to-end style. After alignment, total aligned reads for each gene were counted and saved as the initial expression level. Due to the nonuniform normalization of FPKM value,37 we used a tag-wise dispersion mode with the estimated size factors method38 to normalize the gene transcriptional levels. After normalization, we selected the mRNA and long noncoding RNA (lncRNA) genes whose average transcriptional levels were no less than 0.1 as the expressed genes. 
We then calculated the Z-score for each gene across all samples. To retrieve genes showing diverse expression levels, we set the mean to Z-score 0.7 as the threshold to filter the genes with mediocre expression in all samples. DEGs were identified using the edgeR39 package from R software (R Project for Statistical Computing, Vienna, Austria). The false discovery rate was < 0.05 and |log2 fold change| > 1). 
Correlation Network Analysis
To obtain the expression modules and distinguish genes from a union set by expression features, we used weighted gene coexpression network analysis (WGCNA).40 The normalized expression file of selected genes was the input. The output was the gene modules according to their expression patterns. For each gene module, eigengene was chosen to represent the expression pattern. 
We used Pearson's correlation coefficient (PCC) analysis to investigate the coexpressed genes with VEGFA. We set absolute (PCC) ≥ 0.5 and P values ≤ 0.05 as the thresholds to filter the coexpressed genes. We then obtained the human TFs from AnimalTFDB 2.0,41 and selected the TFs that were coexpressed with VEGFA from the coexpressed genes. 
In Situ Hybridization (ISH) Experiments
ISH was performed using a digoxigenin mRNA detection kit with some modifications (Boster Biological Technology, Wuhan, China). Single-stranded sense and antisense digoxigenin-labeled corresponding mRNA probes were generated by in vitro cDNA transcription with T3 or T7 RNA polymerase. The tissues were fixed in 4% paraformaldehyde (containing 1/1000 DEPC) for 1 hour, then dehydrated in ethanol gradients of 70%, 80%, 90%, 95%, and 100%, and incubated with xylene and embedded in paraffin. Next, 4-mm-thick paraffin sections were baked face-up for 2 hours at 50°C and deparaffinized. The tissue sections were then treated with proteinase K (5 mg/mL; Sigma-Aldrich Corp., St. Louis, MO, USA) in 10 mM Tris-HCl (pH 8.0), 1 mM EDTA at 37°C for 15 minutes before treatment with 3% H2O2 for 8 minutes, then rinsed with phosphate-buffered saline and postfixed in 1% paraformaldehyde containing 1/1000 DEPC at room temperature for 10 minutes. Sections were then prehybridized for 4 hours at 40°C with 50 μL prehybridization buffer. After the prehybridization, hybridization was performed in buffer containing prehybridization fluid and 1.0 mg·L−1 either antisense or sense digoxigenin-labeled mRNA probes at 40°C overnight. Next, the sections were washed at 37°C for 30 minutes in a buffer containing 50% formamide and 2× SSC, followed by digestion with ribonuclease A (10 mg/mL; Sigma-Aldrich Corp.) at 37°C for 30 minutes. After being washed in 2× SSC at 37°C for 20 minutes, and twice in 0.2× SSC at 37°C for 30 minutes, blocking was performed in blocking buffer containing 1% bovine serum albumin for each section for 30 minutes at 37°C. After blocking buffer was removed, the sections were incubated with biotinylated mouse anti-digoxigenin antibody in a humidified chamber at 37°C for 60 minutes, then postrinsed and incubated with Strept Avidin Biotin-peroxidase Complex (SABC) solution for 30 minutes, and finally reacted with a solution containing biotinylated peroxidase for 30 minutes at room temperature. After the sections were rinsed, a coloring reaction was performed with a DAB Color Development Kit (Solarbio Science & Technology, Beijing, China). When the color reaction was visible, it was stopped by washing in buffer solution, and the sections were then counterstained with hematoxylin. Finally, the section was washed, dehydrated, and mounted. 
Western Blotting
Total protein was extracted by cell lysis. The extracts were electroblotted onto a polyvinylidene difluoride (PVDF) membrane following separation by 12% SDS-PAGE. The PVDF membrane was incubated with a blocking solution of 5% skim milk for 1.5 hours at room temperature, followed by incubation with a specific antibody (VEGF-C, 22601-1-AP; Proteintech, Wuhan, China) overnight at 4°C and then with horseradish peroxidase–conjugated secondary antibody (SA00001-2; Proteintech) for 1 hour at room temperature. Finally, the blots were visualized using an ECL reagent (Amersham, Buckinghamshire, UK). 
Statistical Analysis
To describe the biologically functional significance of selected gene sets, we used the hypergeometric test to calculate the enrichment of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Other statistical analysis was performed using R software (R Project for Statistical Computing) unless otherwise stated. The significance of differences was evaluated with either Student's t-test, when only two groups were compared, or the hypergeometric test for a Venn diagram. No statistical methods were used to predetermine sample size. *P < 0.05; **P < 0.01; ***P < 0.001. Hierarchical clustering was performed by Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm, in the public domain) or the heat map function in R. 
Results
Transcriptome Profiles Distinguish NVMs From Retinal Membranes
In order to establish a link between the gene expression profiles in NVMs and PDR, we collected pathological tissues from retinal NVMs of patients that comprehensively included the major types of retinal neovascular disorders. Qualified transcriptome data were obtained from 23 NVM samples from patients, including 17 from T2D patients, 3 from T1D patients, and 3 from BRVO patients (Supplementary Table S1). T2D NVM samples observed at seven different regions were collected (Fig. 1A; Supplementary Table S1). This sampling roughly reflected the incidence of PDR in our hospital. Qualified transcriptome data were also obtained from three retinal membranes from normal individuals and two from T2D patients. A total of 1.43 billion reads were sequenced with an average of 47.51 million reads per sample. After adaptor and low-quality read filtering, approximately 40 million reads per sample were kept and aligned to the human genome (Supplementary Table S2). Similarly to a previous study,42 the sequencing results contained a large fraction of reads (19%–50%) within intronic regions (Supplementary Fig. S1A), and approximately 80% of the reads were not spliced (Supplementary Fig. S1B), indicating that a large number of nascent transcripts were captured. We obtained 25,920 expressed mRNA and lncRNA genes with normalized expressions ≥ 0.1 at least in one sample. 
Figure 1
 
Transcriptome profiles of PDR samples were separated from nonproliferative samples. (A) A fundus photograph shows the right eye of a patient with 20/200 visual acuity, no signs of hypertension. Black arrows show the proliferative membranes (pale yellow) around the vascular arcade. The white arrow shows the blood vessel arcade. Other parts of the retina were labeled on the photograph that was related to Supplementary Table S1. (B) Principal component analysis of the sequenced samples revealed the separation between proliferative membrane and other samples. (C) An UpSet plot showing the overlapped differentially expressed genes between proliferative membrane and normal retina. Only groups with more than 10 genes are shown. DEGs from Ishikawa et al.34 were also included and are shown as a brown color.
Figure 1
 
Transcriptome profiles of PDR samples were separated from nonproliferative samples. (A) A fundus photograph shows the right eye of a patient with 20/200 visual acuity, no signs of hypertension. Black arrows show the proliferative membranes (pale yellow) around the vascular arcade. The white arrow shows the blood vessel arcade. Other parts of the retina were labeled on the photograph that was related to Supplementary Table S1. (B) Principal component analysis of the sequenced samples revealed the separation between proliferative membrane and other samples. (C) An UpSet plot showing the overlapped differentially expressed genes between proliferative membrane and normal retina. Only groups with more than 10 genes are shown. DEGs from Ishikawa et al.34 were also included and are shown as a brown color.
To obtain a global understanding of transcriptional dynamics from NVMs and retinas, we performed principal component analysis (PCA) based on the normalized expression level of all samples (Fig. 1B). NVMs and retinas could be clearly separated by the second principal component (PC2), indicating the presence of a large number of DEGs between NVMs and retinas. NVMs originating from different regions or pathological states could not be clearly separated because of the large individual variation reflected by the PC1 direction. We found that six NVMs (P_I_W, two P_II_VAs, P_II_W, P_BRVO, and P_II_S) were separated from other samples (Fig. 1B). 
We also calculated PCCs for all samples. From the hierarchical clustering results, samples were clustered into three major groups (Supplementary Fig. S1C, blue frames). Group I (top left) contained 11 NVM samples and one retina sample; group II (middle) contained two normal and two T2D retinal samples; and group III (bottom right) contained 12 NVM samples. Group II could be treated as the retina group. Six proliferative samples (P_BRVO3, P_II_S2, P_II_VA1, P_II_W4, P_II_VA2, and P_I_W3) in group I showed low coefficients (R < 0.75) with most of the samples in group III (Supplementary Fig. S1C), consistent with the PCA results (Fig. 1B). The PCA and PCC analyses both showed that the transcriptional profile of NVMs was different from the retinal membranes, consistent with the fact that the NVM was essentially different from the retina.4,43 
We then compared the DEGs between each type of NVM and normal retinas. DEGs in NVMs from the surrounding area of the eye (P_II_S) were more distinct from other groups, and those in the vascular arch and optic disc (P_II_VA_OD) displayed the largest DEG number (3322, Fig. 1C), among which 911 (27.4%) were unique to this group and the rest of DEGs were shared by at least one other group. A total of 1387 DEGs (32.35%) were shared by at least four groups, and 1601 DEGs (37.35%) were identified only from one specific group. A total of 133 DEGs were common in all NVMs versus retina comparisons (Fig. 1C). Only 19 of the 87 DEGs from Ishikawa et al.34 were not included in our DEG set (Fig. 1C), implying the high sensitivity of the RNA-Seq method used in this study. 
Expression of the VEGF Family and Two Other Vascular Endothelium-Specific Growth Factors in NVMs and Retinas
It has been reported that several proteins are major contributors to angiogenesis and vascular development, including five proteins from the VEGF family, four proteins from the angiopoietin family, and at least one protein from the large ephrin family.15,16,44,45 By presenting the transcriptional levels of all these growth factors and their receptors, we found a significant decline of VEGFA and its major receptor KDR in the NVM group (Fig. 2A). This expression pattern suggested that the VEGFA pathway was repressed in NVM samples. Expression of VEGFB did not show significant alterations between NVMs and diabetic retinas. In contrast, VEGFC and its receptor gene, VEGFR-3 (FLT4), exhibited an increased expression pattern (Fig. 2A, left). Both angiopoietin genes and their receptors showed a consistently higher expression level in the NVMs than the diabetic retinas (Fig. 2A, middle). For the ephrin family, expression of EFNB2, EFNA1, and the receptor EPHB4 was increased in the NVMs (Fig. 2A, right). 
Figure 2
 
Expression levels of three gene families were mainly elevated in neovascular membranes. (A) Box plot showing the gene expression level changes between the NVMs and the retinas. Three family genes of vascular growth factors (top) and their receptors (bottom) were presented. The brown box represents the NVM groups and the green box represents the retinal groups. (B) In situ hybridization validation for the expression and location of VEGF family genes. Upper: presentation for P_II_N sample; bottom: for P_II_VA_OD sample. Arrows in the figures represent the locations of target probe signal (dark brown regions). Control: treatment without target probes. (C) Western blotting results showing the protein level of VEGFC in different kinds of NVMs from patients.
Figure 2
 
Expression levels of three gene families were mainly elevated in neovascular membranes. (A) Box plot showing the gene expression level changes between the NVMs and the retinas. Three family genes of vascular growth factors (top) and their receptors (bottom) were presented. The brown box represents the NVM groups and the green box represents the retinal groups. (B) In situ hybridization validation for the expression and location of VEGF family genes. Upper: presentation for P_II_N sample; bottom: for P_II_VA_OD sample. Arrows in the figures represent the locations of target probe signal (dark brown regions). Control: treatment without target probes. (C) Western blotting results showing the protein level of VEGFC in different kinds of NVMs from patients.
To further validate the above observations, we included normal retina in the comparison group. We analyzed the relative expression patterns of some known neovascularization-related genes by sample-based normalized expression levels (Supplementary Fig. S2A). Notably, we found that VEGFA, FLT1, KDR, and tumor neovascularization-related TFs HIF1A and PROX146,47 were highly expressed in normal retinas. In the diabetic state, expression of VEGFA increased dramatically in two T2D retinal samples, while expression of PROX1 was exclusively shut down. Remarkably, VEGFA expression in all NVMs from T2D patients was drastically decreased compared with the T2D retinas (Supplementary Fig. S2A). In contrast, we found that VEGFC was at a low level in normal and T2D retinas, but at a high level in most of the NVM samples (Supplementary Fig. S2A). The differentially regulated expression of VEGFA and VEGFC was also plotted by their gene-based normalized expression levels (Supplementary Fig. S2B). To validate the induced expression of VEGFC, we performed RNA ISH in NVM samples (Fig. 2B). VEGF expression (dark brown) emerged at the margin or the inner membrane (arrows, Fig. 2B) of the tissues. The VEGFC signal was quite strong, similar to that of VEGFA (Fig. 2B, upper). We also validated the increased protein levels of VEGFC in T2D NVMs (Fig. 2C). These results together showed that expression of VEGFC was highly induced in the NVMs, while VEGFA was the most expressed in diabetic retinas. 
Genes Coexpressed With VEGFA and Their Network in the NVMs of PDR Patients
We then characterized the VEGFA-related regulatory network by using a coexpression network approach. A correlation matrix between VEGFA and the other 25,919 genes was generated by computing their PCCs. With criteria of P value ≤ 0.01 and absolute PCCs ≥ 0.5, a total of 1615 genes were coexpressed with VEGFA, including 1421 positive and 194 negative pairs (Supplementary Table S3). Among these coexpressed genes, we found HMGA1, a TF that controls insulin receptor (INSR) gene transcription,48 was positively coexpressed with VEGFA (Fig. 3A), coincident with the previous finding that HMGA1 positively regulated the expression level of VEGFA.49 ANGPT2 was negatively coexpressed with VEGFA (R = −0.729; Fig. 3B), consistent with the upregulation of ANGPT2 in NVMs. In a similar manner as the VEGFA expression patterns, genes positively coexpressed with VEGFA were higher in diabetic and normal retinas than in NVMs (Supplementary Fig. S3A). In contrast, 194 genes that negatively coexpressed with VEGFA showed a reciprocal expression pattern (Supplementary Fig. S3A). 
Figure 3
 
Functions of genes coexpressed with VEGFA and their expression patterns. (A, B) The dot plot represents the coexpression relationship between HMGA1 (A), ANGPT2 (B), and VEGFA. (C) The bar plot of enriched biological process terms by genes positively coexpressed with VEGFA. (D) The bar plot of enriched biological process terms by genes negatively coexpressed with VEGFA. (E) A heat map presentation of the genes belonging to the VEGF signaling pathway. Two group gene sets of the VEGF signaling pathway were labeled by blue and red colors.
Figure 3
 
Functions of genes coexpressed with VEGFA and their expression patterns. (A, B) The dot plot represents the coexpression relationship between HMGA1 (A), ANGPT2 (B), and VEGFA. (C) The bar plot of enriched biological process terms by genes positively coexpressed with VEGFA. (D) The bar plot of enriched biological process terms by genes negatively coexpressed with VEGFA. (E) A heat map presentation of the genes belonging to the VEGF signaling pathway. Two group gene sets of the VEGF signaling pathway were labeled by blue and red colors.
Based on the concept that coexpressed genes are functionally related,50 we further explored the VEGFA-regulated networks by functional enrichment analysis. The genes positively coexpressed with VEGFA were mostly enriched in the glycolysis and glucose metabolic processes (Fig. 3C; Supplementary Fig. S3B). Consistent with the role of hypoxia in increasing VEGFA expression,51 the HIF-1 signaling pathway (KEGG) was also enriched (Supplementary Fig. S3B). The biological process terms, including blood coagulation and platelet activation, as well as regulation of the apoptotic process, and vascular smooth muscle contraction were enriched by genes negatively coexpressed with VEGFA (Fig. 3D; Supplementary Fig. S3C). We also plotted the expression profiles of genes in VEGF signaling pathways including 19 genes. Expression of these genes was classified into two subgroups, with one group (GI) showing a constitutively higher expression in retina but lower expression in NVMs (Fig. 3E). Genes of another group (GII) were expressed at lower levels in normal retina and higher in T2D retina, but expressed at a much higher level in NVMs (Fig. 3E), suggesting that modulation of the VEGF signaling expression network could be a key event for initiating neovascularization. 
A Specific Group of ETS Transcription Factors Promote the Progression of PDR
To further explore VEGFA-mediated gene expression networks in diabetic retinas and NVMs, we analyzed the TFs involved in VEGFA regulation. We found that ETS1 was negatively coexpressed with VEGFA. ETS1 has been recently shown to amplify the VEGF signal through acetylation to promote angiogenesis.52 Two ETSs (ETS1 and ETS2) and two ERGs (ERG and FLI1) domain-containing genes showed a reverse expression pattern with VEGFA, that is, expressing at high levels in NVMs (Fig. 4A). Three of these four genes (ETS1, ERG, and FLI1) are known to promote angiogenesis by activating VEGFR and TIE1/TIE2.53,54 We also found that VEGFA was negatively coexpressed with TIE1 (R = −0.40) and TIE2 (also known as TEK; R = −0.45). By adding VEGFC to this network, we observed the positive coexpression between VEGFC and ETS family members (Fig. 4A). 
Figure 4
 
ETS TFs promote the progression of neovascularization and PDR. (A) Correlation presentation of the four ETS TF genes with VEGFA and VEGFC. VEGFA was negatively coexpressed with these four genes, but VEGFC was positively coexpressed. (B) Hierarchical clustering heat map of the ETS TF family genes. Two distinct clusters emerged. (C) A Venn diagram shows the coexpression genes with the four ETS TF genes. (D) A bar plot of the enriched biological process terms for genes coexpressed with ETS1. (E) A Venn diagram showing the overlapped genes that were coexpressed with ETS1 and bound by ETS1. (F) Regulatory schema plot shows the coexpression relationships between ETS family genes and vascular growth family genes. The interactions in this network were generated by Reactom FI application from Cytoscape (https://cytoscape.org, in the public domain).
Figure 4
 
ETS TFs promote the progression of neovascularization and PDR. (A) Correlation presentation of the four ETS TF genes with VEGFA and VEGFC. VEGFA was negatively coexpressed with these four genes, but VEGFC was positively coexpressed. (B) Hierarchical clustering heat map of the ETS TF family genes. Two distinct clusters emerged. (C) A Venn diagram shows the coexpression genes with the four ETS TF genes. (D) A bar plot of the enriched biological process terms for genes coexpressed with ETS1. (E) A Venn diagram showing the overlapped genes that were coexpressed with ETS1 and bound by ETS1. (F) Regulatory schema plot shows the coexpression relationships between ETS family genes and vascular growth family genes. The interactions in this network were generated by Reactom FI application from Cytoscape (https://cytoscape.org, in the public domain).
It is known that the ETS family contains both transcriptional activators and repressors for angiogenesis.54 To define the global function of ETS family genes in PDR, we searched for all ETS family genes in humans and studied their expression patterns. A clustering heat map of 28 expressed ETS family genes exhibited two distinct expression clusters (Fig. 4B). Cluster I contained 16 genes that showed sporadic expression patterns. Among cluster I genes, ELF2 showed a similar overexpression pattern as VEGFA in the T2D retinas. In contrast, cluster II contained 12 genes, which showed a low expression level in retina samples but globally overexpressed levels in NVM samples (Fig. 4B), implying their potential role in promoting PDR-related neovascularization.54,55 
We then obtained a coexpressed gene network of ETS1/2, EGR, and FLI1 to further explore their potential functions in NVMs. The coexpressed genes with these four TFs significantly overlapped. For example, 732 genes were coexpressed with at least three of these four TFs (26.65%, P = 2.30e-155; hypergeometric test; Fig. 4C). Consistently, functional analysis of ETS-coexpressed genes showed significant enrichment of angiogenesis and related biological processes (Fig. 4D; Supplementary Figs. S4A–C). To assess whether these coexpressed genes resulted from transcriptional activation by TF binding, we downloaded and analyzed the DNA binding profile data of ETS1 from human umbilical vein endothelial cells.52 We found 305 coexpressed genes with ETS1, which were also bound by ETS1 (P = 1.24e-5, hypergeometric test; Fig. 4E), suggesting that a significant fraction of the coexpressed genes with these TFs involved direct transcriptional regulation. 
We then characterized the coexpression relationship between these four TFs and vascular growth factors in Figure 2A. Most notably, besides VEGFC, ANGPT2, EFNB2, and EFNA1, their receptors TIE1/TEK, VEGFR1/3, and EPHB4 were positively coexpressed with all these four TFs (Supplementary Table S4). This implied that the ETS and ERG domains containing TFs could participate in PDR progression by promoting angiogenesis. Based on the coexpression results, we propose that ETS TF family genes, especially ETS1/2, ERG, and FLI1, participated in the formation and progression of PDR and angiogenesis via positively regulating the expression of TIE, VEGFR, and ephrin family genes (Fig. 4F). 
Gene Expression in T2D Retina Is Reprogrammed
We identified gene sets associated with specific sample groups using the WGCNA method. We first calculated the mean Z-score for each gene (see Materials and Methods and Supplementary Fig. S5A) to select genes with a mean Z-score no less than 0.7 (12,293 genes), and then clustered these genes using the WGCNA method40 to obtain a consensus coexpression network. We identified 12 modules with characteristic expression patterns (Supplementary Fig. S5B). These modules could be further classified into two major groups based on clustering of eigengene values (Supplementary Fig. S5C): retinal group (four modules in the right branch), and NVM group (six modules in the left branch). The four retinal modules had higher expression levels in R_norm and R_II samples than others (Figs. 5A–D). Genes in the blue (2341 genes) module were specifically expressed in normal and T2D retinas, and classified as the constitutive retina genes (Fig. 5A). Green (662 genes) and green–yellow modules (99 genes) had higher and specific expression levels in normal retina, which we named T2D-repressed modules (Figs. 5B, 5D). The purple module (214 genes) exhibited consistently high levels in T2D retina samples, which we named T2D-induced genes (Fig. 5C). In summary, the T2D-repressed and -induced genes comprised 29.4% of the four retina groups, indicating that the diabetic state profoundly altered the transcriptional profile of the retina. 
Figure 5
 
WGCNA analysis of the varied genes revealed distinct expression pattern modules. (AD) A bar plot showing the eigengene value (left) and the top 10 enriched biological processes (right) of retinal related modules. Samples were sorted by sample ID order. (E) A heat map presentation showing the top enriched biological process terms for genes in the retina-elated modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
Figure 5
 
WGCNA analysis of the varied genes revealed distinct expression pattern modules. (AD) A bar plot showing the eigengene value (left) and the top 10 enriched biological processes (right) of retinal related modules. Samples were sorted by sample ID order. (E) A heat map presentation showing the top enriched biological process terms for genes in the retina-elated modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
We further explored the functions of genes in the four retina groups using GO and KEGG enrichment analysis. Among the top 10 enriched GO biological process (BP) terms in the blue module, seven of them were directly retinal related, including phototransduction, retinal bipolar neuron differentiation, retina development in camera-type eye, and rhodopsin-mediated signaling pathway (Fig. 5A). This enrichment pattern was highly consistent with the constitutive expression of these genes in the retina, regardless of the pathological state. For enrichment of KEGG pathways, NOD-like receptor, phosphatidylinositol, and neurotrophin, as well as Fc epsilon RI signaling pathways were among the top 10 enriched pathways (Supplementary Fig. S5D). 
T2D-repressed genes were mainly located in the green module. GO enrichment analysis revealed that these genes were highly enriched in synaptic transmission and neuronal functions, including visual perception, sensory perception of sound, and adult walking behavior, as well as response to stimulus (Fig. 5B). KEGG analysis revealed that a large number of synaptic functions were selectively repressed under the T2D state, including serotonergic synapse, cholinergic synapse, GABAergic synapse, and glutamatergic synapse (Supplementary Fig. S5E). In contrast to the T2D-repressed genes, we noticed that the T2D-induced genes in the purple module were specifically enriched in glucose metabolic processes and small molecule metabolic processes, and responses to stress and protein targeting to mitochondrion (Fig. 5C; Supplementary Fig. S5F). T2D-induced genes including VEGFA, ALDOA, ENO1, EGLN1, and LDHA were enriched in the HIF-1 signaling pathway. These results suggested that the transcriptional profiles in T2D retina had adapted to a stressful living condition containing high levels of glucose. 
To compare the differential expression of genes in the retina modules and NVM modules, the top 25 most enriched BP terms and KEGG pathways in the four retina-specific modules were selected and plotted (Fig. 5E; Supplementary Fig. S5G), revealing that T2D-induced retina functions were generally shared more by NVM groups than the T2D repressed. 
The Neovascularization-Specific Gene Expression Programs
We then analyzed NVM-related modules, including eight modules exhibiting high levels in some or most NVM samples, which were called NVM-related modules (7837 genes). The turquoise (4306 genes) and brown (1073 genes) modules exhibited high expression in less than half of the NVMs (Supplementary Figs. S6A–B). The other six modules (2459 genes), including red, black, yellow, pink, magenta and tan, were prevalently high in NVM samples, showing similar or complementary expression patterns (Figs. 6A–F). For example, the red module showed elevated levels of gene expression in most of P_I_W and P_II_W samples and some other NVM samples (Fig. 6A). These genes were named neovascularization-related genes, which could serve as potential PDR markers for further study. We provided detailed module information in the supplementary information section (Supplementary Table S5). 
Figure 6
 
Neovascularization-specific gene expression modules and their enriched functional pathways. (AF) A bar plot showing the eigengene value of six epigene modules showing prevalent upregulated expression in neovascularization samples. Samples were sorted by sample ID order. (G) A heat map presentation of the top enriched biological process terms for genes in each module, choosing terms enriched in neovascularization-specific modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
Figure 6
 
Neovascularization-specific gene expression modules and their enriched functional pathways. (AF) A bar plot showing the eigengene value of six epigene modules showing prevalent upregulated expression in neovascularization samples. Samples were sorted by sample ID order. (G) A heat map presentation of the top enriched biological process terms for genes in each module, choosing terms enriched in neovascularization-specific modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
We then plotted the top 25 BP and KEGG pathways in the six neovascularization-specific modules (Fig. 6G; Supplementary Fig. S6C). Most of these enriched terms were commonly from multiple neovascularization-specific modules. Impressively, most NVM terms were shared among neovascularization-specific modules, but were excluded from the retina groups. The pink module included VEGFC, EPHB4, EFNB1, and other angiogenesis-related genes such as PDGFA, TGFBR1, FZD6, ITGAV, ITGA5, TGFBI, FGFR1, and SHC1. Some angiogenesis growth factors and receptors described in Figure 2 were included in the magenta module, including ANGPT2, EFNA1, TIE2, FLT4, FLT1, and EPHB4 (Supplementary Table S5). Furthermore, the pink (281 genes) and magenta (188 genes) modules were highly enriched in angiogenesis and extracellular matrix organization pathways (Fig. 6G). The expression profile of genes in angiogenesis and the related functional pathways had a global higher expression level in NVM samples (Supplementary Fig. S6D). In addition, several angiogenesis-related pathways such as blood coagulation,56,57 platelet activation,58 and leukocyte migration59 were also highly enriched in NVM modules (Fig. 6G; Supplementary Fig. S6C). 
Discussion
Although PDR is one of the leading causes of preventable blindness, a better understanding of its pathogenetic mechanisms is needed to improve clinical treatments.8 In the present study, we have presented a new approach to study retinal angiogenesis and neovascularization associated with PDR, through robustly comparing the deep-sequenced transcriptomes of NVMs from T2D, T1D, and BRVO patients with those from T2D and normal retinas. To the best of our knowledge, we are the first to characterize the global alterations of transcriptome profiles between NVMs and retinas. The results showed a large transcriptional difference between NVMs and the retina, which provides an unprecedented opportunity to understand the mechanisms of retinal angiogenesis and neovascularization, as well as to develop novel biomarkers and therapies against PDR. 
Gene Expression in T2D Retinas Is Extensively Reprogrammed to Adapt to High-Glucose Stress
It has been reported that loss of neuroretinal adaptation was detected before observable vascular pathologies,60 suggesting preceding transcriptome changes before DR diagnosis. Using the WGCNA method, we identified three modules of retina-specific genes: constitutive retina genes (2341, 70.6%), T2D-repressed retina genes (761, 22.9%), and T2D-induced retina genes (214, 6.5%). Although we obtained only five retinal samples in this study, we could clearly separate the retinal samples from the NVMs to obtain three distinct retina-specific modules. Functional clustering analysis showed that genes in these modules were enriched in either expected or interesting functional pathways, which supported the validity of the results. We observed that the glucose metabolic process and small molecule metabolic process pathways were enriched only in T2D-induced modules. The different pathways between T2D-induced and T2D-repressed modules suggested that the transcriptional profile of retinal cells had been reprogrammed to adapt to a stressful condition including high levels of glucose, and some retinal functions were compromised in T2D patients. These findings are worthy of further investigation in animal models and in the early clinical diagnosis of DR.61 
Genes Overexpressed in NVMs as Potential Biomarkers
Approximately half of patients fail to respond to current anti-VEGF treatment,8 necessitating the study of other factors contributing to PDR to serve as more comprehensive prognostic and diagnostic biomarkers. The study of additional genes contributing to PDR progression is hindered by technical difficulties in studying the very small amounts of patient tissues. A microarray analysis of gene expression in six NVMs from patients with PDR and three normal human retinas has only revealed 91 genes preferentially expressed in active neovascular patients, when compared with normal retinas, excluding VEGFs and angiopoietins.34 In this study, we reported the detection of 1387 overlapped DEGs among four diabetic proliferative membrane groups and the normal retina, including reported active neovascular genes. WGCNA analysis identified six modules with prevalent higher expression in NVM samples, when compared with normal retinas, containing 2459 genes. These genes were highly enriched in angiogenesis, blood coagulation, platelet activation, leukocyte migration, and extracellular matrix organization, consistent with the current understanding of the mechanism of angiogenesis,17,32 and suggested that a large number of angiogenesis-related genes were highly expressed in NVMs. Genes in NVM modules therefore represented potential biomarkers and candidates for developing PDR diagnostic and therapeutic methodologies.62 
Expression of VEGFC, ANGPT, EFNB3, and Their Receptors Is Induced During Retina Neovascularization
The three major families, including VEGFs and the angiopoietin and ephrin protein families, play substantially differential and cooperative roles in promoting and regulating angiogenesis and neovascularization.16,17 In this study, analysis of the expression patterns of all genes in these three families showed that expression of VEGFC, ANGPT1 and ANGPT2, EFNB3, and their receptor genes were coordinately increased in NVMs compared with those of retinas, and these observations were confirmed by ISH and Western blot experiments. The essential role of VEGFC in angiogenesis has been reported in embryonically developing mice and during cancer.6366 The high expression of VEGFC in T2D-induced NVMs suggested its potential role in maintaining the progression of NVM development. These findings represent the first report of a network of three major families of vascular endothelium-specific growth factors that may cooperate during retinal neovascularization, which is consistent with recent reports.27,31 
VEGFA Expression Is Highly Regulated During Retina Pathogenesis and Neovascularization
The essential function of VEGFA in promoting neovascularization in ischemic retinopathy has been proposed based on multiple lines of evidence.9,25,33,67 Our results showed that VEGFA was expressed at high levels in normal retinas, diabetic retinas, and NVMs, while levels of VEGFA in diabetic retinas were significantly higher than in normal retinas and NVMs. The finding that VEGFA expression is specifically increased in T2D retina and reduced in NVMs was unexpected. However, it is consistent with the finding that the onset of DR may be present at the time of diagnoses for T2D patients, and suggests a model in which VEGFA plays an important role for the pathological transformation of retina from T2D to PDR.61 It also provides mechanistic insight into clinical observations that earlier anti-VEGF therapy produces better results.68 
Based on the results presented in this study, we propose a dynamic gene regulatory model that initiates and stabilizes retinal neovascularization and neovascularization. In this model, expression of VEGFA is high in diabetic retina and is sensitive to pathological states. Upon pathological induction, VEGFA secretion into ocular fluid is increased. The increased level of extracellular VEGFA stimulates the expression of ETS TFs in endothelial cells.52 Based on these coexpression relationships, ETS TFs may in turn regulate the expression of growth factors, VEGFC, angiopoietin, and ephrin, as well as their receptors54 (Figs. 2, 6E). These growth factors then coordinately drive the gene expression in endothelial cells to a pattern toward the NVMs rather than retina cells. It is very possible that the relatively high basal levels of VEGFA, and its main receptor, KDR, maintain the proper response of retinal endothelial cells to extracellular VEGFA, and this newly established expression pattern collectively supports the initiation and stabilization of retinal neovascularization. 
In summary, the transcriptome sequencing and bioinformatics analysis presented in this study have not only generated the most informative transcriptome data, a comprehensive list of candidate genes for developing PDR treatment, and biomarkers, but have also provided mechanistic insights into the mechanism of retina angiogenesis and neovascularization. The dynamic expression of VEGFA, the induced expression of endothelium-specific growth factors, and a group of ETS family of TFs form the framework of an elaborate network that coordinates angiogenesis in PDR. Based on these observations, further studies with higher-resolution methods in animal and cell models are needed to further clarify, validate, and extend these regulatory mechanisms that occur during retinal neovascularization in PDR.69 
Acknowledgments
The authors thank Jian Zhou for the discussion about data processing and Hong Wu, ABLife, Inc., for language polishing. 
Supported by the National Natural Science Foundation of China (No. 81272242 and No. 81502284), National Key New Drug Creation and Manufacturing Program of Ministry of Science and Technology (No. 2013ZX09103003004), the Fundamental Research Funds for the Central Universities and the Grant of Jilin Province Science & Technology Committee, China (No. 20170414028GH, No. 20170520031JH, No. 201706230TC, No. 20150204038YY, No. 20150101188JC, No. 20150204038YY, and No. 20150309003YY). This work was also supported by Appreciate the Beauty of Life, Inc. (No. 2014-11032, YZ). 
Disclosure: Y. Li, None; D. Chen, None; L. Sun, None; Y. Wu, None; Y. Zou, None; C. Liang, None; Y. Bao, None; J. Yi, None; Y. Zhang, None; J. Hou, None; Z. Li, None; F. Yu, None; Y. Huang, None; C. Yu, None; L. Liu, None; Z. Liu, None; Y. Zhang, None; Y. Li, None 
References
American Diabetes Association. Diagnosis and classification of diabetes mellitus. Diabetes Care. 2009; 33 (suppl 1): S62–S69.
Saaddine JB, Honeycutt AA, Narayan KM, Zhang X, Klein R, Boyle JP. Projection of diabetic retinopathy and other major eye diseases among people with diabetes mellitus: United States, 2005-2050. Arch Ophthalmol. 2008; 126: 1740–1747.
Giacco F, Brownlee M. Oxidative stress and diabetic complications. Circ Res. 2010; 107: 1058–1070.
Antonetti DA, Klein R, Gardner TW. Diabetic retinopathy. N Engl J Med. 2012; 366: 1227–1239.
Cheung N, Mitchell P, Wong TY. Diabetic retinopathy. Lancet. 2010; 376: 124–136.
Tan GS, Cheung N, Simó R, Cheung GCM, Wong TY. Diabetic macular oedema. Lancet Diabetes Endocrinol. 2017; 5: 143–155.
Lang GE. Diabetic macular edema. Ophthalmologica. 2012; 227 (suppl 1): 21–29.
Stitt AW, Curtis TM, Chen M, et al. The progress in understanding and treatment of diabetic retinopathy. Prog Retin Eye Res. 2016; 51: 156–186.
Ferrara N, Adamis AP. Ten years of anti-vascular endothelial growth factor therapy. Nat Rev Drug Discov. 2016; 15: 385–403.
Leung DW, Cachianes G, Kuang WJ, Goeddel DV, Ferrara N. Vascular endothelial growth factor is a secreted angiogenic mitogen. Science. 1989; 246: 1306–1309.
Senger DR, Galli SJ, Dvorak AM, Perruzzi CA, Harvey VS, Dvorak HF. Tumor cells secrete a vascular permeability factor that promotes accumulation of ascites fluid. Science. 1983; 219: 983–985.
Neufeld G, Cohen T, Gengrinovitch S, Poltorak Z. Vascular endothelial growth factor (VEGF) and its receptors. FASEB J. 1999; 13: 9–22.
Shweiki D, Itin A, Soffer D, Keshet E. Vascular endothelial growth factor induced by hypoxia may mediate hypoxia-initiated angiogenesis. Nature. 1992; 359: 843–845.
Ferrara N, Gerber HP, LeCouter J. The biology of VEGF and its receptors. Nat Med. 2003; 9: 669–676.
Shibuya M. Vascular endothelial growth factor and its receptor system: physiological functions in angiogenesis and pathological roles in various diseases. J Biochem. 2013; 153: 13–19.
Yancopoulos GD, Davis S, Gale NW, Rudge JS, Wiegand SJ, Holash J. Vascular-specific growth factors and blood vessel formation. Nature. 2000; 407: 242–248.
Carmeliet P. Angiogenesis in health and disease. Nat Med. 2003; 9: 653–660.
Lu M, Kuroki M, Amano S, et al. Advanced glycation end products increase retinal vascular endothelial growth factor expression. J Clin Invest. 1998; 101: 1219–1224.
Lu M, Amano S, Miyamoto K, et al. Insulin-induced vascular endothelial growth factor expression in retina. Invest Ophthalmol Vis Sci. 1999; 40: 3281–3286.
Ferrara N. Vascular endothelial growth factor: basic science and clinical progress. Endocr Rev. 2004; 25: 581–611.
Miller JW, Adamis AP, Shima DT, et al. Vascular endothelial growth factor/vascular permeability factor is temporally and spatially correlated with ocular angiogenesis in a primate model. Am J Pathol. 1994; 145: 574–584.
Aiello LP, Avery RL, Arrigg PG, et al. Vascular endothelial growth factor in ocular fluid of patients with diabetic retinopathy and other retinal disorders. N Engl J Med. 1994; 331: 1480–1487.
Adamis AP, Shima DT, Yeo KT, et al. Synthesis and secretion of vascular permeability factor/vascular endothelial growth factor by human retinal pigment epithelial cells. Biochem Biophys Res Commun. 1993; 193: 631–638.
Boulton M, Foreman D, Williams G, Mcleod D. VEGF localisation in diabetic retinopathy. Br J Ophthalmol. 1998; 82: 561–568.
Stewart MW. A review of ranibizumab for the treatment of diabetic retinopathy. Ophthalmol Ther. 2017; 6: 33–47.
Cheung N, Wong IY, Wong TY. Ocular anti-VEGF therapy for diabetic retinopathy: overview of clinical efficacy and evolving applications. Diabetes Care. 2014; 37: 900–905.
Joukov V, Pajusola K, Kaipainen A, et al. A novel vascular endothelial growth factor, VEGF-C, is a ligand for the Flt4 (VEGFR-3) and KDR (VEGFR-2) receptor tyrosine kinases. EMBO J. 1996; 15: 290–298.
Partanen TA, Arola J, Saaristo A, et al. VEGF-C and VEGF-D expression in neuroendocrine cells and their receptor, VEGFR-3, in fenestrated blood vessels in human tissues. FASEB J. 2000; 14: 2087–2096.
Tammela T, Zarkada G, Wallgard E, et al. Blocking VEGFR-3 suppresses angiogenic sprouting and vascular network formation. Nature. 2008; 454: 656–660.
Zarkada G, Heinolainen K, Makinen T, Kubota Y, Alitalo K. VEGFR3 does not sustain retinal angiogenesis without VEGFR2. Proc Natl Acad Sci U S A. 2015; 112: 761–766.
Zhao B, Ma A, Cai J, Boulton M. VEGF-A regulates the expression of VEGF-C in human retinal pigment epithelial cells. Br J Ophthalmol. 2006; 90: 1052–1059.
Carmeliet P. Mechanisms of angiogenesis and arteriogenesis. Nat Med. 2000; 6: 389–395.
Apte RS, Chen DS, Ferrara N. VEGF in signaling and disease: beyond discovery and development. Cell. 2019; 176: 1248–1264.
Ishikawa K, Yoshida S, Kobayashi Y, et al. Microarray analysis of gene expression in fibrovascular membranes excised from patients with proliferative diabetic retinopathy. Invest Ophthalmol Vis Sci. 2015; 56: 932–946.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011; 17: 10–12.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013; 14: R36.
Dillies MA, Rau A, Aubert J, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2013; 14: 671–683.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010; 11: R106.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010; 26: 139–140.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9: 559.
Zhang HM, Liu T, Liu CJ, et al. AnimalTFDB 2.0: a resource for expression, prediction and functional study of animal transcription factors. Nucleic Acids Res. 2015; 43: D76–D81.
Ameur A, Zaghlool A, Halvardson J, et al. Total RNA sequencing reveals nascent transcription and widespread co-transcriptional splicing in the human brain. Nat Struct Mol Biol. 2011; 18: 1435–1440.
Ljubimov AV, Burgeson RE, Butkowski RJ, et al. Basement membrane abnormalities in human eyes with diabetic retinopathy. J Histochem Cytochem. 1996; 44: 1469–1479.
Loukovaara S, Gucciardo E, Repo P, et al. Indications of lymphatic endothelial differentiation and endothelial progenitor cell activation in the pathology of proliferative diabetic retinopathy. Acta Ophthalmol. 2015; 93: 512–523.
Domigan CK, Ziyad S, Iruela-Arispe ML. Canonical and noncanonical vascular endothelial growth factor pathways: new developments in biology and signal transduction. Arterioscler Thromb Vasc Biol. 2015; 35: 30–39.
Hillen F, Griffioen AW. Tumour vascularization: sprouting angiogenesis and beyond. Cancer Metastasis Rev. 2007; 26: 489–502.
Semenza GL. Targeting HIF-1 for cancer therapy. Nat Rev Cancer. 2003; 3: 721–732.
Foti D, Iuliano R, Chiefari E, Brunetti A. A nucleoprotein complex containing Sp1, C/EBP beta, and HMGI-Y controls human insulin receptor gene transcription. Mol Cell Biol. 2003; 23: 2720–2732.
Chiefari E, Ventura V, Capula C, et al. A polymorphism of HMGA1 protects against proliferative diabetic retinopathy by impairing HMGA1-induced VEGFA expression. Sci Rep. 2016; 6: 39429.
Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci U S A. 1998; 95: 14863–14868.
Simo R, Sundstrom JM, Antonetti DA. Ocular anti-VEGF therapy for diabetic retinopathy: the role of VEGF in the pathogenesis of diabetic retinopathy. Diabetes Care. 2014; 37: 893–899.
Chen J, Fu Y, Day DS, et al. VEGF amplifies transcription through ETS1 acetylation to enable angiogenesis. Nat Commun. 2017; 8: 383.
Sharrocks AD. The ETS-domain transcription factor family. Nat Rev Mol Cell Biol. 2001; 2: 827–837.
Lelievre E, Lionneton F, Soncin F, Vandenbunder B. The Ets family contains transcriptional activators and repressors involved in angiogenesis. Int J Biochem Cell Biol. 2001; 33: 391–407.
Watanabe D, Takagi H, Suzuma K, et al. Transcription factor Ets-1 mediates ischemia- and vascular endothelial growth factor-dependent retinal neovascularization. Am J Pathol. 2004; 164: 1827–1835.
Belting M, Ahamed J, Ruf W. Signaling of the tissue factor coagulation pathway in angiogenesis and cancer. Arterioscler Thromb Vasc Biol. 2005; 25: 1545–1550.
Rothmeier AS, Liu E, Chakrabarty S, et al. Identification of the coagulation factor VIIa integrin binding site required for pro-angiogenic PAR2 signaling. Blood. 2018; 131: 674–685.
Janowska-Wieczorek A, Wysoczynski M, Kijowski J, et al. Microvesicles derived from activated platelets induce metastasis and angiogenesis in lung cancer. Int J Cancer. 2005; 113: 752–760.
Kreuger J, Phillipson M. Targeting vascular and leukocyte communication in angiogenesis, inflammation and fibrosis. Nat Rev Drug Discov. 2016; 15: 125–142.
Abcouwer SF, Gardner TW. Diabetic retinopathy: loss of neuroretinal adaptation to the diabetic metabolic environment. Ann N Y Acad Sci. 2014; 1311: 174–190.
Ebneter A, Zinkernagel MS. Novelties in diabetic retinopathy. Endocr Dev. 2016; 31: 84–96.
MacLaren RE, Bennett J, Schwartz SD. Gene therapy and stem cell transplantation in retinal disease: the new frontier. Ophthalmology. 2016; 123: S98–S106.
Bauer SM, Bauer RJ, Liu ZJ, Chen H, Goldstein L, Velazquez OC. Vascular endothelial growth factor-C promotes vasculogenesis, angiogenesis, and collagen constriction in three-dimensional collagen gels. J Vasc Surg. 2005; 41: 699–707.
Hamada K, Oike Y, Takakura N, et al. VEGF-C signaling pathways through VEGFR-2 and VEGFR-3 in vasculoangiogenesis and hematopoiesis. Blood. 2000; 96: 3793–3800.
Bauer SM, Bauer RJ, Velazquez OC. Angiogenesis, vasculogenesis, and induction of healing in chronic wounds. Vasc Endovasc Surg. 2005; 39: 293–306.
Valtola R, Salven P, Heikkila P, et al. VEGFR-3 and its ligand VEGF-C are associated with angiogenesis in breast cancer. Am J Pathol. 1999; 154: 1381–1390.
Miller JW. VEGF: from discovery to therapy: the Champalimaud Award Lecture. Trans Vis Sci Tech. 2016; 5 (2): 9.
Brown DM, Nguyen QD, Marcus DM, et al. Long-term outcomes of ranibizumab therapy for diabetic macular edema: the 36-month results from two phase III trials: RISE and RIDE. Ophthalmology. 2013; 120: 2013–2022.
Gariano RF, Gardner TW. Retinal angiogenesis in development and disease. Nature. 2005; 438: 960–966.
Figure 1
 
Transcriptome profiles of PDR samples were separated from nonproliferative samples. (A) A fundus photograph shows the right eye of a patient with 20/200 visual acuity, no signs of hypertension. Black arrows show the proliferative membranes (pale yellow) around the vascular arcade. The white arrow shows the blood vessel arcade. Other parts of the retina were labeled on the photograph that was related to Supplementary Table S1. (B) Principal component analysis of the sequenced samples revealed the separation between proliferative membrane and other samples. (C) An UpSet plot showing the overlapped differentially expressed genes between proliferative membrane and normal retina. Only groups with more than 10 genes are shown. DEGs from Ishikawa et al.34 were also included and are shown as a brown color.
Figure 1
 
Transcriptome profiles of PDR samples were separated from nonproliferative samples. (A) A fundus photograph shows the right eye of a patient with 20/200 visual acuity, no signs of hypertension. Black arrows show the proliferative membranes (pale yellow) around the vascular arcade. The white arrow shows the blood vessel arcade. Other parts of the retina were labeled on the photograph that was related to Supplementary Table S1. (B) Principal component analysis of the sequenced samples revealed the separation between proliferative membrane and other samples. (C) An UpSet plot showing the overlapped differentially expressed genes between proliferative membrane and normal retina. Only groups with more than 10 genes are shown. DEGs from Ishikawa et al.34 were also included and are shown as a brown color.
Figure 2
 
Expression levels of three gene families were mainly elevated in neovascular membranes. (A) Box plot showing the gene expression level changes between the NVMs and the retinas. Three family genes of vascular growth factors (top) and their receptors (bottom) were presented. The brown box represents the NVM groups and the green box represents the retinal groups. (B) In situ hybridization validation for the expression and location of VEGF family genes. Upper: presentation for P_II_N sample; bottom: for P_II_VA_OD sample. Arrows in the figures represent the locations of target probe signal (dark brown regions). Control: treatment without target probes. (C) Western blotting results showing the protein level of VEGFC in different kinds of NVMs from patients.
Figure 2
 
Expression levels of three gene families were mainly elevated in neovascular membranes. (A) Box plot showing the gene expression level changes between the NVMs and the retinas. Three family genes of vascular growth factors (top) and their receptors (bottom) were presented. The brown box represents the NVM groups and the green box represents the retinal groups. (B) In situ hybridization validation for the expression and location of VEGF family genes. Upper: presentation for P_II_N sample; bottom: for P_II_VA_OD sample. Arrows in the figures represent the locations of target probe signal (dark brown regions). Control: treatment without target probes. (C) Western blotting results showing the protein level of VEGFC in different kinds of NVMs from patients.
Figure 3
 
Functions of genes coexpressed with VEGFA and their expression patterns. (A, B) The dot plot represents the coexpression relationship between HMGA1 (A), ANGPT2 (B), and VEGFA. (C) The bar plot of enriched biological process terms by genes positively coexpressed with VEGFA. (D) The bar plot of enriched biological process terms by genes negatively coexpressed with VEGFA. (E) A heat map presentation of the genes belonging to the VEGF signaling pathway. Two group gene sets of the VEGF signaling pathway were labeled by blue and red colors.
Figure 3
 
Functions of genes coexpressed with VEGFA and their expression patterns. (A, B) The dot plot represents the coexpression relationship between HMGA1 (A), ANGPT2 (B), and VEGFA. (C) The bar plot of enriched biological process terms by genes positively coexpressed with VEGFA. (D) The bar plot of enriched biological process terms by genes negatively coexpressed with VEGFA. (E) A heat map presentation of the genes belonging to the VEGF signaling pathway. Two group gene sets of the VEGF signaling pathway were labeled by blue and red colors.
Figure 4
 
ETS TFs promote the progression of neovascularization and PDR. (A) Correlation presentation of the four ETS TF genes with VEGFA and VEGFC. VEGFA was negatively coexpressed with these four genes, but VEGFC was positively coexpressed. (B) Hierarchical clustering heat map of the ETS TF family genes. Two distinct clusters emerged. (C) A Venn diagram shows the coexpression genes with the four ETS TF genes. (D) A bar plot of the enriched biological process terms for genes coexpressed with ETS1. (E) A Venn diagram showing the overlapped genes that were coexpressed with ETS1 and bound by ETS1. (F) Regulatory schema plot shows the coexpression relationships between ETS family genes and vascular growth family genes. The interactions in this network were generated by Reactom FI application from Cytoscape (https://cytoscape.org, in the public domain).
Figure 4
 
ETS TFs promote the progression of neovascularization and PDR. (A) Correlation presentation of the four ETS TF genes with VEGFA and VEGFC. VEGFA was negatively coexpressed with these four genes, but VEGFC was positively coexpressed. (B) Hierarchical clustering heat map of the ETS TF family genes. Two distinct clusters emerged. (C) A Venn diagram shows the coexpression genes with the four ETS TF genes. (D) A bar plot of the enriched biological process terms for genes coexpressed with ETS1. (E) A Venn diagram showing the overlapped genes that were coexpressed with ETS1 and bound by ETS1. (F) Regulatory schema plot shows the coexpression relationships between ETS family genes and vascular growth family genes. The interactions in this network were generated by Reactom FI application from Cytoscape (https://cytoscape.org, in the public domain).
Figure 5
 
WGCNA analysis of the varied genes revealed distinct expression pattern modules. (AD) A bar plot showing the eigengene value (left) and the top 10 enriched biological processes (right) of retinal related modules. Samples were sorted by sample ID order. (E) A heat map presentation showing the top enriched biological process terms for genes in the retina-elated modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
Figure 5
 
WGCNA analysis of the varied genes revealed distinct expression pattern modules. (AD) A bar plot showing the eigengene value (left) and the top 10 enriched biological processes (right) of retinal related modules. Samples were sorted by sample ID order. (E) A heat map presentation showing the top enriched biological process terms for genes in the retina-elated modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
Figure 6
 
Neovascularization-specific gene expression modules and their enriched functional pathways. (AF) A bar plot showing the eigengene value of six epigene modules showing prevalent upregulated expression in neovascularization samples. Samples were sorted by sample ID order. (G) A heat map presentation of the top enriched biological process terms for genes in each module, choosing terms enriched in neovascularization-specific modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
Figure 6
 
Neovascularization-specific gene expression modules and their enriched functional pathways. (AF) A bar plot showing the eigengene value of six epigene modules showing prevalent upregulated expression in neovascularization samples. Samples were sorted by sample ID order. (G) A heat map presentation of the top enriched biological process terms for genes in each module, choosing terms enriched in neovascularization-specific modules. Enrichment scores (−log10 P value) were normalized by terms for the presentation.
Supplement 1
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×